Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 2 Mar 2018 23:10:08 +0000 (23:10 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 2 Mar 2018 23:10:08 +0000 (23:10 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15397 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 97031e8..9a25348 100644 (file)
@@ -149,7 +149,7 @@ public abstract class CVODES {
         * ----------------------------------------\r
         */\r
 \r
-       final int  CV_SUCCESS = 0;\r
+       final int CV_SUCCESS = 0;\r
        final int CV_TSTOP_RETURN = 1;\r
        final int CV_ROOT_RETURN = 2;\r
 \r
@@ -217,6 +217,7 @@ public abstract class CVODES {
        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 HUB_FACTOR = 0.1;\r
        final double HLB_FACTOR = 100.0;\r
@@ -231,11 +232,37 @@ public abstract class CVODES {
     final double RDIV = TWO;\r
     // MSBP  max no. of steps between lsetup calls\r
     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
+    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
+    //             reset eta =  SUNMIN(eta, ETAMXF)\r
+    final double MXNEF1 = 3;\r
+    final double SMALL_NEF = 2;\r
+\r
+    final double ETACF = 0.25;\r
+    final double ADDON  = 0.000001;\r
+    final double BIAS1 = 6.0;\r
+    final double BIAS2 = 6.0;\r
+    final double BIAS3 = 10.0;\r
+    // LONG_WAIT   number of steps to wait before considering an order change when\r
+    //             q==1 and MXNEF1 error test failures have occurred\r
+    final int LONG_WAIT = 10;\r
+\r
 \r
-       \r
        final int RTFOUND = +1;\r
        final int CLOSERT = +3;\r
 \r
+       /*\r
+        * Control constants for sensitivity DQ\r
+        * ------------------------------------\r
+        */\r
+\r
+       final int CENTERED1 = +1;\r
+       final int CENTERED2 = +2;\r
+       final int FORWARD1 = +3;\r
+       final int FORWARD2 = +4;\r
 \r
        \r
        final String MSG_TIME = "t = %lg";\r
@@ -758,22 +785,9 @@ 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
        \r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
@@ -987,8 +1001,9 @@ public abstract class CVODES {
          double cv_acnrmQS;         /* | acorQS |                                  */\r
          double cv_nlscoef;         /* coeficient in nonlinear convergence test    */\r
          int  cv_mnewt;               /* Newton iteration counter                    */\r
-         int  cv_ncfS1[];              /* Array of Ns local counters for conv.  \r
+         int  cv_ncfS1[][];              /* Array of Ns local counters for conv.  \r
                                        * failures (used in CVStep for STAGGERED1)    */\r
+                                     /* int[][1]\r
 \r
          /*------\r
            Limits \r
@@ -1023,16 +1038,17 @@ public abstract class CVODES {
          long cv_nfQeS;       /* number of fQ calls from sensi DQ                */\r
 \r
 \r
-         long cv_ncfn;        /* number of corrector convergence failures        */\r
-         long cv_ncfnS;       /* number of total sensi. corr. conv. failures     */\r
-         long cv_ncfnS1[];     /* number of sensi. corrector conv. failures       */\r
+         long cv_ncfn[] = new long[1];        /* number of corrector convergence failures        */\r
+         long cv_ncfnS[] = new long[1];       /* number of total sensi. corr. conv. failures     */\r
+         long cv_ncfnS1[][];     /* number of sensi. corrector conv. failures       */\r
+                                 // long[][1]\r
 \r
          long cv_nni;         /* number of nonlinear iterations performed        */\r
          long cv_nniS;        /* number of total sensi. nonlinear iterations     */\r
          long cv_nniS1[];      /* number of sensi. nonlinear iterations           */\r
 \r
-         long cv_netf;        /* number of error test failures                   */\r
-         long cv_netfQ;       /* number of quadr. error test failures            */\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
 \r
@@ -1502,8 +1518,8 @@ public abstract class CVODES {
 \r
        cv_mem.cv_nst     = 0;\r
        cv_mem.cv_nfe     = 0;\r
-       cv_mem.cv_ncfn    = 0;\r
-       cv_mem.cv_netf    = 0;\r
+       cv_mem.cv_ncfn[0]    = 0;\r
+       cv_mem.cv_netf[0]    = 0;\r
        cv_mem.cv_nni     = 0;\r
        cv_mem.cv_nsetups = 0;\r
        cv_mem.cv_nhnil   = 0;\r
@@ -2230,7 +2246,7 @@ public abstract class CVODES {
       }\r
 \r
       if (cv_mem.cv_quadr_sensi) {\r
-        retval = fQS(cv_mem.cv_Ns, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+        retval = cvQuadSensRhsInternalDQ(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
@@ -2666,10 +2682,10 @@ public abstract class CVODES {
     /* 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
+      cv_mem.cv_ncfnS[0] = 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
+        cv_mem.cv_ncfnS[0] += cv_mem.cv_ncfnS1[is][0];\r
       }\r
     }\r
     \r
@@ -3711,12 +3727,12 @@ public abstract class CVODES {
    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
+     retval = cvSensRhsInternalDQ(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
+       retval = cvSensRhs1InternalDQ(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
@@ -4046,7 +4062,7 @@ public abstract class CVODES {
    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
+     retval = cvQuadSensRhsInternalDQ(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
@@ -4817,11 +4833,18 @@ else                return(snrm);
 \r
  private int cvStep(CVodeMemRec cv_mem)\r
  {\r
-   double saved_t, dsm, dsmQ, dsmS, dsmQS;\r
+   double dsm[] = new double[1];\r
+   double saved_t;\r
+   double dsmQ[] = new double[1];\r
+   double 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 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 nflag[] = new int [1];\r
+   int kflag, eflag;\r
    int retval, is;\r
 \r
    /* Are we computing sensitivities with a staggered approach? */\r
@@ -4831,12 +4854,13 @@ else                return(snrm);
 \r
    /* Initialize local counters for convergence and error test failures */\r
 \r
-   ncf  = nef  = 0;\r
-   nefQ = nefQS = 0;\r
-   ncfS = nefS = 0;\r
+   ncf[0]  = nef[0]  = 0;\r
+   nefQ[0] = nefQS = 0;\r
+   ncfS[0] = nefS = 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
-       cv_mem.cv_ncfS1[is] = 0;\r
+       cv_mem.cv_ncfS1[is][0] = 0;\r
    }\r
 \r
    /* If needed, adjust method parameters */\r
@@ -4847,7 +4871,7 @@ else                return(snrm);
    /* Looping point for attempts to take a step */\r
 \r
    saved_t = cv_mem.cv_tn;\r
-   nflag = FIRST_CALL;\r
+   nflag[0] = FIRST_CALL;\r
 \r
    for(;;) {  \r
 \r
@@ -4856,92 +4880,92 @@ else                return(snrm);
 \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
+     nflag[0] = cvNls(cv_mem, nflag[0]);\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
+     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
+     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
+     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
+     if (eflag == TRY_AGAIN) continue;\r
 \r
      /* Return if error test failed and recovery not possible. */\r
-    /* if (eflag != CV_SUCCESS) return(eflag); */\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
+     if (cv_mem.cv_quadr) {\r
 \r
-       ncf = nef = 0;*/ /* reset counters for states */\r
+       ncf[0] = nef[0] = 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
+       nflag[0] = 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
+       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
+       if (cv_mem.cv_errconQ) {\r
+         cv_mem.cv_acnrmQ = N_VWrmsNorm_Serial(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
+         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
+         if (dsmQ[0] > dsm[0]) dsm[0] = dsmQ[0];\r
        }\r
 \r
-     }*/\r
+     }\r
 \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 = nef = 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; */ /* 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
+     /* 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
          continue;\r
        }\r
 \r
-       if (do_sensi_stg) {*/\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
+        /* nflag[0] = 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
+         /*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
            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
@@ -6292,129 +6316,1532 @@ else                return(snrm);
 \r
        /* Call the lsolve function */\r
        b = cv_mem.cv_tempv;\r
-       /*retval = cv_mem->cv_lsolve(cv_mem, b, cv_mem.cv_ewt,\r
-                                  cv_mem.cv_y, cv_mem.cv_ftemp); \r
-       cv_mem->cv_nni++;\r
+       retval = cv_lsolve(cv_mem, b, cv_mem.cv_ewt,\r
+                                  cv_mem.cv_y, cv_mem.cv_ftemp, cv_mem.cv_lsolve_select); \r
+       cv_mem.cv_nni++;\r
 \r
-       if (retval < 0) return(CV_LSOLVE_FAIL);*/\r
+       if (retval < 0) return(CV_LSOLVE_FAIL);\r
        \r
        /* If lsolve had a recoverable failure and Jacobian data is\r
           not current, signal to try the solution again            */\r
-      /* if (retval > 0) { \r
-         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+       if (retval > 0) { \r
+         if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
            return(TRY_AGAIN);\r
          else\r
            return(CONV_FAIL);\r
-       }*/\r
+       }\r
 \r
        /* Solve the sensitivity linear systems and do the same \r
           tests on the return value of lsolve. */\r
     \r
-      /* if (do_sensi_sim) {\r
+       if (do_sensi_sim) {\r
 \r
-         for (is=0; is<cv_mem->cv_Ns; is++) {\r
-           N_VLinearSum(cv_mem->cv_rl1, cv_mem->cv_znS[1][is], ONE,\r
-                        cv_mem->cv_acorS[is], cv_mem->cv_tempvS[is]);\r
-           N_VLinearSum(cv_mem->cv_gamma, cv_mem->cv_ftempS[is], -ONE,\r
-                        cv_mem->cv_tempvS[is], cv_mem->cv_tempvS[is]);\r
+         for (is=0; is<cv_mem.cv_Ns; is++) {\r
+           N_VLinearSum_Serial(cv_mem.cv_rl1, cv_mem.cv_znS[1][is], ONE,\r
+                        cv_mem.cv_acorS[is], cv_mem.cv_tempvS[is]);\r
+           N_VLinearSum_Serial(cv_mem.cv_gamma, cv_mem.cv_ftempS[is], -ONE,\r
+                        cv_mem.cv_tempvS[is], cv_mem.cv_tempvS[is]);\r
          }\r
-         bS = cv_mem->cv_tempvS;\r
-         for (is=0; is<cv_mem->cv_Ns; is++) {\r
-           retval = cv_mem->cv_lsolve(cv_mem, bS[is], cv_mem->cv_ewtS[is],\r
-                                      cv_mem->cv_y, cv_mem->cv_ftemp);\r
+         bS = cv_mem.cv_tempvS;\r
+         for (is=0; is<cv_mem.cv_Ns; is++) {\r
+           retval = cv_lsolve(cv_mem, bS[is], cv_mem.cv_ewtS[is],\r
+                                      cv_mem.cv_y, cv_mem.cv_ftemp, cv_mem.cv_lsolve_select);\r
            if (retval < 0) return(CV_LSOLVE_FAIL);\r
            if (retval > 0) { \r
-             if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+             if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
                return(TRY_AGAIN);\r
              else\r
                return(CONV_FAIL);\r
            }\r
          }\r
-       }*/\r
+       }\r
        \r
        /* Get WRMS norm of correction; add correction to acor and y */\r
 \r
-       /*del = N_VWrmsNorm(b, cv_mem->cv_ewt);\r
-       N_VLinearSum(ONE, cv_mem->cv_acor, ONE, b, cv_mem->cv_acor);\r
-       N_VLinearSum(ONE, cv_mem->cv_zn[0], ONE, cv_mem->cv_acor, cv_mem->cv_y);\r
+       del = N_VWrmsNorm_Serial(b, cv_mem.cv_ewt);\r
+       N_VLinearSum_Serial(ONE, cv_mem.cv_acor, ONE, b, cv_mem.cv_acor);\r
+       N_VLinearSum_Serial(ONE, cv_mem.cv_zn[0], ONE, cv_mem.cv_acor, cv_mem.cv_y);\r
 \r
        if (do_sensi_sim) {\r
-         delS = cvSensUpdateNorm(cv_mem, del, bS, cv_mem->cv_ewtS);\r
-         for (is=0; is<cv_mem->cv_Ns; is++) {\r
-           N_VLinearSum(ONE, cv_mem->cv_acorS[is], ONE,\r
-                        bS[is], cv_mem->cv_acorS[is]);\r
-           N_VLinearSum(ONE, cv_mem->cv_znS[0][is], ONE,\r
-                        cv_mem->cv_acorS[is], cv_mem->cv_yS[is]);\r
+         delS = cvSensUpdateNorm(cv_mem, del, bS, cv_mem.cv_ewtS);\r
+         for (is=0; is<cv_mem.cv_Ns; is++) {\r
+           N_VLinearSum_Serial(ONE, cv_mem.cv_acorS[is], ONE,\r
+                        bS[is], cv_mem.cv_acorS[is]);\r
+           N_VLinearSum_Serial(ONE, cv_mem.cv_znS[0][is], ONE,\r
+                        cv_mem.cv_acorS[is], cv_mem.cv_yS[is]);\r
          }\r
-       }*/\r
+       }\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
-       /*Del = (do_sensi_sim) ? delS : del;\r
+       Del = (do_sensi_sim) ? delS : del;\r
        if (m > 0)\r
-         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
+         cv_mem.cv_crate = Math.max(CRDOWN * cv_mem.cv_crate, Del/Delp);\r
+       dcon = Del * Math.min(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
+           if (do_sensi_sim && cv_mem.cv_errconS)\r
+             cv_mem.cv_acnrm = delS;\r
            else\r
-             cv_mem->cv_acnrm = del;\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
+           cv_mem.cv_acnrm = N_VWrmsNorm_Serial(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
-         cv_mem->cv_jcur = SUNFALSE;\r
-         return(CV_SUCCESS); */ /* Convergence achieved */\r
-       /*}\r
+         cv_mem.cv_jcur[0] = false;\r
+         return(CV_SUCCESS);  /* Convergence achieved */\r
+       }\r
 \r
-       cv_mem->cv_mnewt = ++m;*/\r
+       cv_mem.cv_mnewt = ++m;\r
        \r
        /* Stop at maxcor iterations or if iter. seems to be diverging.\r
           If still not converged and Jacobian data is not current, \r
           signal to try the solution again                            */\r
-       /*if ((m == cv_mem->cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp))) {\r
-         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+       if ((m == cv_mem.cv_maxcor) || ((m >= 2) && (Del > RDIV * Delp))) {\r
+         if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
            return(TRY_AGAIN);\r
          else\r
            return(CONV_FAIL);\r
-       }*/\r
+       }\r
        \r
        /* Save norm of correction, evaluate f, and loop again */\r
-      /* Delp = Del;\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
+       Delp = Del;\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
-         if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+         if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
            return(TRY_AGAIN);\r
          else\r
            return(RHSFUNC_RECVR);\r
        }\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_y,\r
-                                   cv_mem->cv_ftemp, cv_mem->cv_yS,\r
-                                   cv_mem->cv_ftempS, wrk1, wrk2);\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_y,\r
+                                   cv_mem.cv_ftemp, cv_mem.cv_yS,\r
+                                   cv_mem.cv_ftempS, wrk1, wrk2);\r
          if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
          if (retval > 0) {\r
-           if ((!cv_mem->cv_jcur) && (cv_mem->cv_lsetup))\r
+           if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
              return(TRY_AGAIN);\r
            else\r
              return(SRHSFUNC_RECVR);\r
          }\r
-       }*/\r
+       }\r
 \r
      } /* end loop */\r
+     \r
+     \r
 \r
    }\r
+   \r
+   private int cv_lsolve(CVodeMemRec cv_mem, NVector b, NVector weight,\r
+                  NVector ycur, NVector fcur, int select) {\r
+          int flag = 0;\r
+          switch (select) {\r
+          case cvDlsSolve_select:\r
+                  flag = cvDlsSolve(cv_mem, b, weight, ycur, fcur);\r
+                  break;\r
+          }\r
+          return flag;\r
+   }\r
+   \r
+   /*-----------------------------------------------------------------\r
+   cvDlsSolve\r
+   -----------------------------------------------------------------\r
+   This routine interfaces between CVode and the generic \r
+   SUNLinearSolver object LS, by calling the solver and scaling \r
+   the solution appropriately when gamrat != 1.\r
+   -----------------------------------------------------------------*/\r
+ private int cvDlsSolve(CVodeMemRec cv_mem, NVector b, NVector weight,\r
+                NVector ycur, NVector fcur)\r
+ {\r
+   int retval;\r
+   CVDlsMemRec cvdls_mem;\r
+\r
+   /* Return immediately if cv_mem or cv_mem->cv_lmem are NULL */\r
+   if (cv_mem == null) {\r
+     cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS", \r
+                   "cvDlsSolve", MSGD_CVMEM_NULL);\r
+     return(CVDLS_MEM_NULL);\r
+   }\r
+   if (cv_mem.cv_lmem == null) {\r
+     cvProcessError(cv_mem, CVDLS_LMEM_NULL, "CVSDLS", \r
+                   "cvDlsSolve", MSGD_LMEM_NULL);\r
+     return(CVDLS_LMEM_NULL);\r
+   }\r
+   cvdls_mem = cv_mem.cv_lmem;\r
+\r
+   /* call the generic linear system solver, and copy b to x */\r
+   retval = SUNLinSolSolve_Dense(cvdls_mem.LS, cvdls_mem.A, cvdls_mem.x, b, ZERO);\r
+   N_VScale_Serial(ONE, cvdls_mem.x, b);\r
+   \r
+   /* scale the correction to account for change in gamma */\r
+   if ((cv_mem.cv_lmm == CV_BDF) && (cv_mem.cv_gamrat != ONE))\r
+     N_VScale_Serial(TWO/(ONE + cv_mem.cv_gamrat), b, b);\r
+   \r
+   /* store solver return value and return */\r
+   cvdls_mem.last_flag = retval;\r
+   return(retval);\r
+ }\r
+\r
+\r
+        private int SUNLinSolSolve_Dense(SUNLinearSolver S, double A[][], NVector x, \r
+                NVector b, double tol)\r
+       {\r
+       double A_cols[][];\r
+    double xdata[];\r
+       int pivots[];\r
+       int i,j;\r
+       \r
+       if ( (A == null) || (S == null) || (x == null) || (b == null) ) \r
+       return(SUNLS_MEM_NULL);\r
+       \r
+       /* copy b into x */\r
+       N_VScale_Serial(ONE, b, x);\r
+       \r
+       /* access data pointers (return with failure on NULL) */\r
+       A_cols = null;\r
+       xdata = null;\r
+       pivots = null;\r
+       A_cols = new double[A[0].length][A.length];\r
+    for (i = 0; i < A.length; i++) {\r
+        for (j = 0; j < A[0].length; j++) {\r
+                A_cols[j][i] = A[i][j];\r
+        }\r
+    }\r
+       xdata = x.data;\r
+       pivots = S.pivots;\r
+       if ( (A_cols == null) || (xdata == null)  || (pivots == null) ) {\r
+       S.last_flag = SUNLS_MEM_FAIL;\r
+       return(S.last_flag);\r
+       }\r
+       \r
+       /* solve using LU factors */\r
+       denseGETRS(A_cols, A.length, pivots, xdata);\r
+       for (i = 0; i < A.length; i++) {\r
+        for (j = 0; j < A[0].length; j++) {\r
+                A[i][j] = A_cols[j][i];\r
+        }\r
+    }\r
+    S.last_flag = SUNLS_SUCCESS;\r
+       return(S.last_flag);\r
+       }\r
+        \r
+        private void denseGETRS(double a[][], int n, int p[], double b[])\r
+        {\r
+          int i, k, pk,q;\r
+          double col_k[] = new double[a[0].length];\r
+          double tmp;\r
+\r
+          /* Permute b, based on pivot information in p */\r
+          for (k=0; k<n; k++) {\r
+            pk = p[k];\r
+            if(pk != k) {\r
+              tmp = b[k];\r
+              b[k] = b[pk];\r
+              b[pk] = tmp;\r
+            }\r
+          }\r
+\r
+          /* Solve Ly = b, store solution y in b */\r
+          for (k=0; k<n-1; k++) {\r
+                for (q = 0; q < a[0].length; q++) {\r
+                col_k[q] = a[k][q];\r
+                }\r
+            for (i=k+1; i<n; i++) b[i] -= col_k[i]*b[k];\r
+          }\r
+\r
+          /* Solve Ux = y, store solution x in b */\r
+          for (k = n-1; k > 0; k--) {\r
+            for (q = 0; q < a[0].length; q++) {\r
+                col_k[q] = a[k][q];\r
+                }\r
+            b[k] /= col_k[k];\r
+            for (i=0; i<k; i++) b[i] -= col_k[i]*b[k];\r
+          }\r
+          b[0] /= a[0][0];\r
+\r
+        }\r
+\r
+        /*\r
+         * cvHandleNFlag\r
+         *\r
+         * This routine takes action on the return value nflag = *nflagPtr\r
+         * returned by cvNls, as follows:\r
+         *\r
+         * If cvNls succeeded in solving the nonlinear system, then\r
+         * cvHandleNFlag returns the constant DO_ERROR_TEST, which tells cvStep\r
+         * to perform the error test.\r
+         *\r
+         * If the nonlinear system was not solved successfully, then ncfn and\r
+         * ncf = *ncfPtr are incremented and Nordsieck array zn is restored.\r
+         *\r
+         * If the solution of the nonlinear system failed due to an\r
+         * unrecoverable failure by setup, we return the value CV_LSETUP_FAIL.\r
+         * \r
+         * If it failed due to an unrecoverable failure in solve, then we return\r
+         * the value CV_LSOLVE_FAIL.\r
+         *\r
+         * If it failed due to an unrecoverable failure in rhs, then we return\r
+         * the value CV_RHSFUNC_FAIL.\r
+         *\r
+         * If it failed due to an unrecoverable failure in quad rhs, then we return\r
+         * the value CV_QRHSFUNC_FAIL.\r
+         *\r
+         * If it failed due to an unrecoverable failure in sensi rhs, then we return\r
+         * the value CV_SRHSFUNC_FAIL.\r
+         *\r
+         * Otherwise, a recoverable failure occurred when solving the \r
+         * nonlinear system (cvNls returned nflag = CONV_FAIL, RHSFUNC_RECVR, or\r
+         * SRHSFUNC_RECVR). \r
+         * In this case, if ncf is now equal to maxncf or |h| = hmin, \r
+         * we return the value CV_CONV_FAILURE (if nflag=CONV_FAIL), or\r
+         * CV_REPTD_RHSFUNC_ERR (if nflag=RHSFUNC_RECVR), or CV_REPTD_SRHSFUNC_ERR\r
+         * (if nflag=SRHSFUNC_RECVR).\r
+         * If not, we set *nflagPtr = PREV_CONV_FAIL and return the value\r
+         * PREDICT_AGAIN, telling cvStep to reattempt the step.\r
+         *\r
+         */\r
+\r
+        private int cvHandleNFlag(CVodeMemRec cv_mem, int nflagPtr[], double saved_t,\r
+                                 int ncfPtr[], long ncfnPtr[])\r
+        {\r
+          int nflag;\r
+          \r
+          nflag = nflagPtr[0];\r
+          \r
+          if (nflag == CV_SUCCESS) return(DO_ERROR_TEST);\r
+\r
+          /* The nonlinear soln. failed; increment ncfn and restore zn */\r
+          (ncfnPtr[0])++;\r
+          cvRestore(cv_mem, saved_t);\r
+          \r
+          /* Return if lsetup, lsolve, or some rhs failed unrecoverably */\r
+          if (nflag == CV_LSETUP_FAIL)    return(CV_LSETUP_FAIL);\r
+          if (nflag == CV_LSOLVE_FAIL)    return(CV_LSOLVE_FAIL);\r
+          if (nflag == CV_RHSFUNC_FAIL)   return(CV_RHSFUNC_FAIL);\r
+          if (nflag == CV_QRHSFUNC_FAIL)  return(CV_QRHSFUNC_FAIL);\r
+          if (nflag == CV_SRHSFUNC_FAIL)  return(CV_SRHSFUNC_FAIL);\r
+          if (nflag == CV_QSRHSFUNC_FAIL) return(CV_QSRHSFUNC_FAIL);\r
+\r
+          /* At this point, nflag = CONV_FAIL, RHSFUNC_RECVR, or SRHSFUNC_RECVR; \r
+             increment ncf */\r
+          \r
+          (ncfPtr[0])++;\r
+          cv_mem.cv_etamax = ONE;\r
+\r
+          /* If we had maxncf failures or |h| = hmin, \r
+             return CV_CONV_FAILURE, CV_REPTD_RHSFUNC_ERR, \r
+             CV_REPTD_QRHSFUNC_ERR, or CV_REPTD_SRHSFUNC_ERR */\r
+\r
+          if ((Math.abs(cv_mem.cv_h) <= cv_mem.cv_hmin*ONEPSM) ||\r
+              (ncfPtr[0] == cv_mem.cv_maxncf)) {\r
+            if (nflag == CONV_FAIL)       return(CV_CONV_FAILURE);\r
+            if (nflag == RHSFUNC_RECVR)   return(CV_REPTD_RHSFUNC_ERR);    \r
+            if (nflag == QRHSFUNC_RECVR)  return(CV_REPTD_QRHSFUNC_ERR);    \r
+            if (nflag == SRHSFUNC_RECVR)  return(CV_REPTD_SRHSFUNC_ERR);    \r
+            if (nflag == QSRHSFUNC_RECVR) return(CV_REPTD_QSRHSFUNC_ERR);    \r
+          }\r
+\r
+          /* Reduce step size; return to reattempt the step */\r
+\r
+          cv_mem.cv_eta = Math.max(ETACF, cv_mem.cv_hmin / Math.abs(cv_mem.cv_h));\r
+          nflagPtr[0] = PREV_CONV_FAIL;\r
+          cvRescale(cv_mem);\r
+\r
+          return(PREDICT_AGAIN);\r
+        }\r
+        \r
+        /*\r
+         * cvRestore\r
+         *\r
+         * This routine restores the value of cv_mem->cv_tn to saved_t and undoes the\r
+         * prediction.  After execution of cvRestore, the Nordsieck array zn has\r
+         * the same values as before the call to cvPredict.\r
+         */\r
+\r
+        private void cvRestore(CVodeMemRec cv_mem, double saved_t)\r
+        {\r
+          int j, k;\r
+          int is;\r
+\r
+          cv_mem.cv_tn = saved_t;\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
+         * Error Test\r
+         * -----------------------------------------------------------------\r
+         */\r
+\r
+        /*\r
+         * cvDoErrorTest\r
+         *\r
+         * This routine performs the local error test, for the state, quadrature, \r
+         * or sensitivity variables. Its last three arguments change depending\r
+         * on which variables the error test is to be performed on.\r
+         * \r
+         * The weighted local error norm dsm is loaded into *dsmPtr, and \r
+         * the test dsm ?<= 1 is made.\r
+         *\r
+         * If the test passes, cvDoErrorTest returns CV_SUCCESS. \r
+         *\r
+         * If the test fails, we undo the step just taken (call cvRestore) and \r
+         *\r
+         *   - if maxnef error test failures have occurred or if SUNRabs(h) = hmin,\r
+         *     we return CV_ERR_FAILURE.\r
+         *\r
+         *   - if more than MXNEF1 error test failures have occurred, an order\r
+         *     reduction is forced. If already at order 1, restart by reloading \r
+         *     zn from scratch (also znQ and znS if appropriate).\r
+         *     If f() fails, we return CV_RHSFUNC_FAIL or CV_UNREC_RHSFUNC_ERR;\r
+         *     if fQ() fails, we return CV_QRHSFUNC_FAIL or CV_UNREC_QRHSFUNC_ERR;\r
+         *     if cvSensRhsWrapper() fails, we return CV_SRHSFUNC_FAIL or CV_UNREC_SRHSFUNC_ERR;\r
+         *     (no recovery is possible at this stage).\r
+         *\r
+         *   - otherwise, set *nflagPtr to PREV_ERR_FAIL, and return TRY_AGAIN. \r
+         *\r
+         */\r
+\r
+        private int cvDoErrorTest(CVodeMemRec cv_mem, int nflagPtr[], double saved_t, \r
+                                 double acor_nrm,\r
+                                 int nefPtr[], long netfPtr[], double dsmPtr[])\r
+        {\r
+          double dsm;\r
+          int retval, is;\r
+          NVector wrk1, wrk2;\r
+\r
+          dsm = acor_nrm * cv_mem.cv_tq[2];\r
+\r
+          /* If est. local error norm dsm passes test, return CV_SUCCESS */  \r
+          dsmPtr[0] = dsm; \r
+          if (dsm <= ONE) return(CV_SUCCESS);\r
+          \r
+          /* Test failed; increment counters, set nflag, and restore zn array */\r
+          (nefPtr[0])++;\r
+          (netfPtr[0])++;\r
+          nflagPtr[0] = PREV_ERR_FAIL;\r
+          cvRestore(cv_mem, saved_t);\r
+\r
+          /* At maxnef failures or |h| = hmin, return CV_ERR_FAILURE */\r
+          if ((Math.abs(cv_mem.cv_h) <= cv_mem.cv_hmin*ONEPSM) ||\r
+              (nefPtr[0] == cv_mem.cv_maxnef))\r
+            return(CV_ERR_FAILURE);\r
+\r
+          /* Set etamax = 1 to prevent step size increase at end of this step */\r
+          cv_mem.cv_etamax = ONE;\r
+\r
+          /* Set h ratio eta from dsm, rescale, and return for retry of step */\r
+          if (nefPtr[0] <= MXNEF1) {\r
+            cv_mem.cv_eta = ONE / (Math.pow(BIAS2*dsm,ONE/cv_mem.cv_L) + ADDON);\r
+            cv_mem.cv_eta = Math.max(ETAMIN, Math.max(cv_mem.cv_eta,\r
+                                                   cv_mem.cv_hmin / Math.abs(cv_mem.cv_h)));\r
+            if (nefPtr[0] >= SMALL_NEF)\r
+              cv_mem.cv_eta = Math.min(cv_mem.cv_eta, ETAMXF);\r
+            cvRescale(cv_mem);\r
+            return(TRY_AGAIN);\r
+          }\r
+          \r
+          /* After MXNEF1 failures, force an order reduction and retry step */\r
+          if (cv_mem.cv_q > 1) {\r
+            cv_mem.cv_eta = Math.max(ETAMIN, cv_mem.cv_hmin / Math.abs(cv_mem.cv_h));\r
+            cvAdjustOrder(cv_mem,-1);\r
+            cv_mem.cv_L = cv_mem.cv_q;\r
+            cv_mem.cv_q--;\r
+            cv_mem.cv_qwait = cv_mem.cv_L;\r
+            cvRescale(cv_mem);\r
+            return(TRY_AGAIN);\r
+          }\r
+\r
+          /* If already at order 1, restart: reload zn, znQ, znS, znQS from scratch */\r
+          cv_mem.cv_eta = Math.max(ETAMIN, cv_mem.cv_hmin / Math.abs(cv_mem.cv_h));\r
+          cv_mem.cv_h *= 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_qwait = LONG_WAIT;\r
+          cv_mem.cv_nscon = 0;\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(CV_UNREC_RHSFUNC_ERR);\r
+\r
+          N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_tempv, cv_mem.cv_zn[1]);\r
+\r
+          if (cv_mem.cv_quadr) {\r
+\r
+            retval = fQ(cv_mem.cv_tn, cv_mem.cv_zn[0],\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(CV_UNREC_QRHSFUNC_ERR);\r
+\r
+            N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_tempvQ, cv_mem.cv_znQ[1]);\r
+\r
+          }\r
+\r
+          if (cv_mem.cv_sensi) {\r
+\r
+            wrk1 = cv_mem.cv_ftemp;\r
+            wrk2 = cv_mem.cv_ftempS[0];\r
+\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(CV_UNREC_SRHSFUNC_ERR);\r
+\r
+            for (is=0; is<cv_mem.cv_Ns; is++) \r
+              N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_tempvS[is], cv_mem.cv_znS[1][is]);\r
+\r
+          }\r
+\r
+          if (cv_mem.cv_quadr_sensi) {\r
+\r
+            wrk1 = cv_mem.cv_ftemp;\r
+            wrk2 = cv_mem.cv_ftempQ;\r
+\r
+            retval = cvQuadSensRhsInternalDQ(cv_mem.cv_Ns, cv_mem.cv_tn,\r
+                                    cv_mem.cv_zn[0], cv_mem.cv_znS[0], \r
+                                    cv_mem.cv_tempvQ, cv_mem.cv_tempvQS,\r
+                                    cv_mem.cv_fQS_data, wrk1, wrk2);\r
+            cv_mem.cv_nfQSe++;\r
+            if (retval < 0) return(CV_QSRHSFUNC_FAIL);\r
+            if (retval > 0) return(CV_UNREC_QSRHSFUNC_ERR);\r
+\r
+            for (is=0; is<cv_mem.cv_Ns; is++) \r
+              N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_tempvQS[is], cv_mem.cv_znQS[1][is]);\r
+\r
+          }\r
+\r
+          \r
+          return(TRY_AGAIN);\r
+        }\r
+\r
+        /*\r
+         * cvQuadNls\r
+         * \r
+         * This routine solves for the quadrature variables at the new step.\r
+         * It does not solve a nonlinear system, but rather updates the\r
+         * quadrature variables. The name for this function is just for \r
+         * uniformity purposes.\r
+         *\r
+         * Possible return values (interpreted by cvHandleNFlag)\r
+         *\r
+         *   CV_SUCCESS       -> continue with error test\r
+         *   CV_QRHSFUNC_FAIL -> halt the integration \r
+         *   QRHSFUNC_RECVR   -> predict again or stop if too many\r
+         *   \r
+         */\r
+\r
+        private int cvQuadNls(CVodeMemRec cv_mem)\r
+        {\r
+          int retval;\r
+\r
+          /* Save quadrature correction in acorQ */\r
+          retval = fQ(cv_mem.cv_tn, cv_mem.cv_y,\r
+                                 cv_mem.cv_acorQ, 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
+          /* If needed, save the value of yQdot = fQ into ftempQ\r
+           * for use in evaluating fQS */\r
+          if (cv_mem.cv_quadr_sensi) {\r
+            N_VScale_Serial(ONE, cv_mem.cv_acorQ, cv_mem.cv_ftempQ);\r
+          }\r
+\r
+          N_VLinearSum_Serial(cv_mem.cv_h, cv_mem.cv_acorQ, -ONE,\r
+                       cv_mem.cv_znQ[1], cv_mem.cv_acorQ);\r
+          N_VScale_Serial(cv_mem.cv_rl1, cv_mem.cv_acorQ, cv_mem.cv_acorQ);\r
+\r
+          /* Apply correction to quadrature variables */\r
+          N_VLinearSum_Serial(ONE, cv_mem.cv_znQ[0], ONE, cv_mem.cv_acorQ, cv_mem.cv_yQ);\r
+\r
+          return(CV_SUCCESS);\r
+        }\r
+        \r
+        /*\r
+         * cvStgrNls\r
+         *\r
+         * This is a high-level routine that attempts to solve the \r
+         * sensitivity linear systems using nonlinear iterations (CV_FUNCTIONAL\r
+         * or CV_NEWTON - depending on the value of iter) once the states y_n\r
+         * were obtained and passed the error test.\r
+         */\r
+\r
+        private int cvStgrNls(CVodeMemRec cv_mem)\r
+        {\r
+          int flag=CV_SUCCESS;\r
+\r
+          switch(cv_mem.cv_iter) {\r
+          case CV_FUNCTIONAL:\r
+            flag = cvStgrNlsFunctional(cv_mem);\r
+            break;\r
+          case CV_NEWTON:\r
+            flag = cvStgrNlsNewton(cv_mem);\r
+            break;\r
+          }\r
+\r
+          return(flag);\r
+\r
+        }\r
+\r
+        /*\r
+         * cvStgrNlsfunctional\r
+         *\r
+         * This routine attempts to solve the sensitivity linear systems \r
+         * using functional iteration (no matrices involved).\r
+         *\r
+         * Possible return values:\r
+         *  CV_SUCCESS,\r
+         *  CV_SRHSFUNC_FAIL,\r
+         *  CONV_FAIL, SRHSFUNC_RECVR\r
+         */\r
+\r
+        private int cvStgrNlsFunctional(CVodeMemRec cv_mem)\r
+        {\r
+          int m;\r
+          int retval, is;\r
+          double Del, Delp, dcon;\r
+          NVector wrk1, wrk2;\r
+\r
+          /* Initialize estimated conv. rate and counter */\r
+          cv_mem.cv_crateS = ONE;\r
+          m = 0;\r
+\r
+          /* Initialize Delp to avoid compiler warning message */\r
+          Delp = ZERO;\r
+\r
+          /* Evaluate fS at predicted yS but with converged y (and corresponding f) */\r
+          wrk1 = cv_mem.cv_tempv;\r
+          wrk2 = cv_mem.cv_ftempS[0];\r
+          retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                    cv_mem.cv_ftemp, 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
+          /* Initialize correction to zero */\r
+          for (is=0; is<cv_mem.cv_Ns; is++)\r
+            N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+\r
+          /* Loop until convergence; accumulate corrections in acorS */\r
+\r
+          for(;;) {\r
+            \r
+            cv_mem.cv_nniS++;\r
+            \r
+            /* Correct yS from last fS value */\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],\r
+                       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
+            /* Get norm of current correction to use in convergence test */\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
+            Del = cvSensNorm(cv_mem, cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
+            for (is=0; is<cv_mem.cv_Ns; is++)\r
+              N_VScale_Serial(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 crateS, and used in the test. \r
+               acnrmS contains the norm of the corrections (yS_n-yS_n(0)) and\r
+               will be used in the error test (if errconS==SUNTRUE)              */\r
+            if (m > 0)\r
+              cv_mem.cv_crateS = Math.max(CRDOWN * cv_mem.cv_crateS, Del / Delp);\r
+            dcon = Del * Math.min(ONE, cv_mem.cv_crateS) / cv_mem.cv_tq[4];\r
+            \r
+            if (dcon <= ONE) {\r
+              if (cv_mem.cv_errconS)\r
+                cv_mem.cv_acnrmS = (m==0)?\r
+                  Del : cvSensNorm(cv_mem, cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
+              return(CV_SUCCESS);  /* Convergence achieved */\r
+            }\r
+\r
+            /* Stop at maxcor iterations or if iter. seems to be diverging */\r
+            m++;\r
+            if ((m==cv_mem.cv_maxcorS) || ((m >= 2) && (Del > RDIV * Delp)))\r
+              return(CONV_FAIL);\r
+\r
+            /* Save norm of correction, evaluate f, and loop again */\r
+            Delp = Del;\r
+\r
+            wrk1 = cv_mem.cv_tempv;\r
+            wrk2 = cv_mem.cv_ftempS[0];\r
+            retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                      cv_mem.cv_ftemp, 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
+          } /* end loop */\r
+          \r
+        }\r
+\r
+        /*\r
+         * cvStgrNlsNewton\r
+         *\r
+         * This routine attempts to solve the sensitivity linear systems using \r
+         * Newton iteration. It calls cvStgrNlsNewton to perform the actual \r
+         * iteration. If the Newton iteration fails with out-of-date Jacobian \r
+         * data (ier=TRY_AGAIN), it calls lsetup and retries the Newton iteration. \r
+         * This second try is unlikely to happen when using a Krylov linear solver.\r
+         *\r
+         * Possible return values:\r
+         *\r
+         *   CV_SUCCESS\r
+         *\r
+         *   CV_LSOLVE_FAIL   -+\r
+         *   CV_LSETUP_FAIL    |\r
+         *   CV_SRHSFUNC_FAIL -+\r
+         *\r
+         *   CONV_FAIL        -+\r
+         *   SRHSFUNC_RECVR   -+\r
+         */\r
+\r
+        private int cvStgrNlsNewton(CVodeMemRec cv_mem)\r
+        {\r
+          int retval, is;\r
+          int convfail, ier;\r
+          NVector vtemp1, vtemp2, vtemp3, wrk1, wrk2;\r
+\r
+          for(;;) {\r
+\r
+            /* Set acorS to zero and load prediction into yS vector */\r
+            for (is=0; is<cv_mem.cv_Ns; is++) {\r
+              N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+              N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], cv_mem.cv_yS[is]);\r
+            }\r
+         \r
+            /* Evaluate fS at predicted yS but with converged y (and corresponding f) */\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_y,\r
+                                      cv_mem.cv_ftemp, cv_mem.cv_yS,\r
+                                      cv_mem.cv_ftempS, wrk1, wrk2);\r
+            if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+            if (retval > 0) return(SRHSFUNC_RECVR);\r
+\r
+            /* Do the Newton iteration */\r
+            ier = cvStgrNewtonIteration(cv_mem);\r
+\r
+            /* If the solve was successful (ier=CV_SUCCESS) or if an error \r
+               that cannot be fixed by a call to lsetup occured\r
+               (ier = CV_LSOLVE_FAIL or CONV_FAIL) return */\r
+            if (ier != TRY_AGAIN) return(ier);\r
+\r
+            /* There was a convergence failure and the Jacobian-related data\r
+               appears not to be current. Call lsetup with convfail=CV_FAIL_BAD_J\r
+               and then loop again */\r
+            convfail = CV_FAIL_BAD_J;\r
+\r
+            /* Rename some vectors for readibility */\r
+            vtemp1 = cv_mem.cv_tempv;\r
+            vtemp2 = cv_mem.cv_yS[0];\r
+            vtemp3 = cv_mem.cv_ftempS[0];\r
+\r
+            /* Call linear solver setup at converged y */\r
+            ier = cv_lsetup(cv_mem, convfail, cv_mem.cv_y,\r
+                                    cv_mem.cv_ftemp, cv_mem.cv_jcur, \r
+                                    vtemp1, vtemp2, vtemp3, cv_mem.cv_lsetup_select);\r
+            cv_mem.cv_nsetups++;\r
+            cv_mem.cv_nsetupsS++;\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
+\r
+            /* Return if lsetup failed */\r
+            if (ier < 0) return(CV_LSETUP_FAIL);\r
+            if (ier > 0) return(CONV_FAIL);\r
+\r
+          } /* end loop */\r
+\r
+        }\r
+\r
+        /*\r
+         * cvStgrNewtonIteration\r
+         *\r
+         * This routine performs the Newton iteration for all sensitivities. \r
+         * If the iteration succeeds, it returns the value CV_SUCCESS. \r
+         * If not, it may signal the cvStgrNlsNewton routine to call lsetup and \r
+         * reattempt the iteration, by returning the value TRY_AGAIN. (In this case, \r
+         * cvStgrNlsNewton must set convfail to CV_FAIL_BAD_J before calling setup again). \r
+         * Otherwise, this routine returns one of the appropriate values \r
+         * CV_LSOLVE_FAIL or CONV_FAIL back to cvStgrNlsNewton.\r
+         */\r
+\r
+        private int cvStgrNewtonIteration(CVodeMemRec cv_mem)\r
+        {\r
+          int m, retval;\r
+          double Del, Delp, dcon;\r
+          NVector bS[];\r
+          NVector wrk1, wrk2;\r
+          int is;\r
+\r
+          m = 0;\r
+\r
+          /* Initialize Delp to avoid compiler warning message */\r
+          Delp = ZERO;\r
+\r
+          /* ftemp  <- f(t_n, y_n)\r
+             y      <- y_n\r
+             ftempS <- fS(t_n, y_n(0), s_n(0))\r
+             acorS  <- 0\r
+             yS     <- yS_n(0)                   */\r
+\r
+          for(;;) {\r
+\r
+            /* Evaluate the residual of the nonlinear systems */\r
+            for (is=0; is<cv_mem.cv_Ns; is++) {\r
+              N_VLinearSum_Serial(cv_mem.cv_rl1, cv_mem.cv_znS[1][is], ONE,\r
+                           cv_mem.cv_acorS[is], cv_mem.cv_tempvS[is]);\r
+              N_VLinearSum_Serial(cv_mem.cv_gamma, cv_mem.cv_ftempS[is], -ONE,\r
+                           cv_mem.cv_tempvS[is], cv_mem.cv_tempvS[is]);\r
+            }\r
+\r
+            /* Call the lsolve function */\r
+            bS = cv_mem.cv_tempvS;\r
+            cv_mem.cv_nniS++;\r
+            for (is=0; is<cv_mem.cv_Ns; is++) {\r
+\r
+              retval = cv_lsolve(cv_mem, bS[is], cv_mem.cv_ewtS[is],\r
+                                         cv_mem.cv_y, cv_mem.cv_ftemp,cv_mem.cv_lsolve_select);\r
+\r
+              /* Unrecoverable error in lsolve */\r
+              if (retval < 0) return(CV_LSOLVE_FAIL);\r
+\r
+              /* Recoverable error in lsolve and Jacobian data not current */\r
+              if (retval > 0) { \r
+                if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                  return(TRY_AGAIN);\r
+                else\r
+                  return(CONV_FAIL);\r
+              }\r
+\r
+            }\r
+         \r
+            /* Get norm of correction; add correction to acorS and yS */\r
+            Del = cvSensNorm(cv_mem, bS, cv_mem.cv_ewtS);\r
+            for (is=0; is<cv_mem.cv_Ns; is++) {\r
+              N_VLinearSum_Serial(ONE, cv_mem.cv_acorS[is], ONE,\r
+                           bS[is], cv_mem.cv_acorS[is]);\r
+              N_VLinearSum_Serial(ONE, cv_mem.cv_znS[0][is], ONE,\r
+                           cv_mem.cv_acorS[is], cv_mem.cv_yS[is]);\r
+            }\r
+\r
+            /* Test for convergence.  If m > 0, an estimate of the convergence\r
+               rate constant is stored in crateS, and used in the test.        */\r
+            if (m > 0)\r
+              cv_mem.cv_crateS = Math.max(CRDOWN * cv_mem.cv_crateS, Del/Delp);\r
+            dcon = Del * Math.min(ONE, cv_mem.cv_crateS) / cv_mem.cv_tq[4];\r
+            if (dcon <= ONE) {\r
+              if (cv_mem.cv_errconS)\r
+                cv_mem.cv_acnrmS = (m==0) ?\r
+                  Del : cvSensNorm(cv_mem, cv_mem.cv_acorS, cv_mem.cv_ewtS);\r
+              cv_mem.cv_jcur[0] = false;\r
+              return(CV_SUCCESS);  /* Convergence achieved */\r
+            }\r
+\r
+            m++;\r
+\r
+            /* Stop at maxcor iterations or if iter. seems to be diverging.\r
+               If still not converged and Jacobian data is not current, \r
+               signal to try the solution again                            */\r
+            if ((m == cv_mem.cv_maxcorS) || ((m >= 2) && (Del > RDIV * Delp))) {\r
+              if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                return(TRY_AGAIN);\r
+              else\r
+                return(CONV_FAIL);\r
+            }\r
+            \r
+            /* Save norm of correction, evaluate fS, and loop again */\r
+            Delp = Del;\r
+            \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_y,\r
+                                      cv_mem.cv_ftemp, cv_mem.cv_yS,\r
+                                      cv_mem.cv_ftempS, wrk1, wrk2);\r
+            if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+            if (retval > 0) {\r
+              if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                return(TRY_AGAIN);\r
+              else\r
+                return(SRHSFUNC_RECVR);\r
+            }\r
+\r
+          } /* end loop */\r
+\r
+        }\r
+\r
+        /*\r
+         * cvStgr1Nls\r
+         *\r
+         * This is a high-level routine that attempts to solve the i-th \r
+         * sensitivity linear system using nonlinear iterations (CV_FUNCTIONAL\r
+         * or CV_NEWTON - depending on the value of iter) once the states y_n\r
+         * were obtained and passed the error test.\r
+         */\r
+\r
+        private int cvStgr1Nls(CVodeMemRec cv_mem, int is)\r
+        {\r
+          int flag=CV_SUCCESS;\r
+\r
+          switch(cv_mem.cv_iter) {\r
+          case CV_FUNCTIONAL:\r
+            flag = cvStgr1NlsFunctional(cv_mem,is);\r
+            break;\r
+          case CV_NEWTON:\r
+            flag = cvStgr1NlsNewton(cv_mem,is);\r
+            break;\r
+          }\r
+\r
+          return(flag);\r
+\r
+        }\r
+\r
+        /*\r
+         * cvStgr1NlsFunctional\r
+         *\r
+         * This routine attempts to solve the i-th sensitivity linear system\r
+         * using functional iteration (no matrices involved).\r
+         *\r
+         * Possible return values:\r
+         *   CV_SUCCESS,\r
+         *   CV_SRHSFUNC_FAIL,\r
+         *   CONV_FAIL, SRHSFUNC_RECVR\r
+         */\r
+\r
+        private int cvStgr1NlsFunctional(CVodeMemRec cv_mem, int is)\r
+        {\r
+          int retval, m;\r
+          double Del, Delp, dcon;\r
+          NVector wrk1, wrk2;\r
+\r
+          /* Initialize estimated conv. rate and counter */\r
+          cv_mem.cv_crateS = ONE;\r
+          m = 0;\r
+\r
+          /* Initialize Delp to avoid compiler warning message */\r
+          Delp = ZERO;\r
+\r
+          /* Evaluate fS at predicted yS but with converged y (and corresponding f) */\r
+          wrk1 = cv_mem.cv_tempv;\r
+          wrk2 = cv_mem.cv_ftempS[0];\r
+          retval = cvSensRhs1Wrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                     cv_mem.cv_ftemp, is, cv_mem.cv_znS[0][is],\r
+                                     cv_mem.cv_tempvS[is], wrk1, wrk2);\r
+          if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+          if (retval > 0) return(SRHSFUNC_RECVR);\r
+\r
+          /* Initialize correction to zero */\r
+          N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+\r
+          /* Loop until convergence; accumulate corrections in acorS */\r
+\r
+          for(;;) {\r
+\r
+            cv_mem.cv_nniS1[is]++;\r
+\r
+            /* Correct yS from last fS value */\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
+            /* Get WRMS norm of current correction to use in convergence test */\r
+            N_VLinearSum_Serial(ONE, cv_mem.cv_tempvS[is], -ONE,\r
+                         cv_mem.cv_acorS[is], cv_mem.cv_acorS[is]);\r
+            Del = N_VWrmsNorm_Serial(cv_mem.cv_acorS[is], cv_mem.cv_ewtS[is]);\r
+            N_VScale_Serial(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 crateS, and used in the test. */\r
+\r
+            if (m > 0)\r
+              cv_mem.cv_crateS = Math.max(CRDOWN * cv_mem.cv_crateS, Del / Delp);\r
+            dcon = Del * Math.min(ONE, cv_mem.cv_crateS) / cv_mem.cv_tq[4];\r
+\r
+            if (dcon <= ONE) {\r
+              return(CV_SUCCESS);  /* Convergence achieved */\r
+            }\r
+\r
+            /* Stop at maxcor iterations or if iter. seems to be diverging */\r
+            m++;\r
+            if ((m==cv_mem.cv_maxcorS) || ((m >= 2) && (Del > RDIV * Delp)))\r
+              return(CONV_FAIL);\r
+\r
+            /* Save norm of correction, evaluate f, and loop again */\r
+            Delp = Del;\r
+\r
+            wrk1 = cv_mem.cv_tempv;\r
+            wrk2 = cv_mem.cv_ftempS[0];\r
+\r
+            retval = cvSensRhs1Wrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                       cv_mem.cv_ftemp, is, cv_mem.cv_yS[is], \r
+                                       cv_mem.cv_tempvS[is], wrk1, wrk2);\r
+            if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+            if (retval > 0) return(SRHSFUNC_RECVR);\r
+\r
+          } /* end loop */\r
+          \r
+        }\r
+\r
+        /*\r
+         * cvStgr1NlsNewton\r
+         *\r
+         * This routine attempts to solve the i-th sensitivity linear system \r
+         * using Newton iteration. It calls cvStgr1NlsNewton to perform the \r
+         * actual iteration. If the Newton iteration fails with out-of-date \r
+         * Jacobian data (ier=TRY_AGAIN), it calls lsetup and retries the \r
+         * Newton iteration. This second try is unlikely to happen when \r
+         * using a Krylov linear solver.\r
+         *\r
+         * Possible return values:\r
+         *\r
+         *   CV_SUCCESS\r
+         *\r
+         *   CV_LSOLVE_FAIL\r
+         *   CV_LSETUP_FAIL\r
+         *   CV_SRHSFUNC_FAIL\r
+         *\r
+         *   CONV_FAIL\r
+         *   SRHSFUNC_RECVR\r
+         */\r
+\r
+        private int cvStgr1NlsNewton(CVodeMemRec cv_mem, int is)\r
+        {\r
+          int convfail, retval, ier;\r
+          NVector vtemp1, vtemp2, vtemp3, wrk1, wrk2;\r
+          \r
+          for(;;) {\r
+\r
+            /* Set acorS to zero and load prediction into yS vector */\r
+            N_VConst_Serial(ZERO, cv_mem.cv_acorS[is]);\r
+            N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], cv_mem.cv_yS[is]);\r
+         \r
+            /* Evaluate fS at predicted yS but with converged y (and corresponding f) */\r
+            wrk1 = cv_mem.cv_tempv;\r
+            wrk2 = cv_mem.cv_tempvS[0];\r
+            retval = cvSensRhs1Wrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                       cv_mem.cv_ftemp, is, cv_mem.cv_yS[is],\r
+                                       cv_mem.cv_ftempS[is], wrk1, wrk2);\r
+            if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+            if (retval > 0) return(SRHSFUNC_RECVR);\r
+          \r
+            /* Do the Newton iteration */\r
+            ier = cvStgr1NewtonIteration(cv_mem, is);\r
+\r
+            /* If the solve was successful (ier=CV_SUCCESS) or if an error \r
+               that cannot be fixed by a call to lsetup occured\r
+               (ier = CV_LSOLVE_FAIL or CONV_FAIL) return */\r
+            if (ier != TRY_AGAIN) return(ier);\r
+\r
+            /* There was a convergence failure and the Jacobian-related data\r
+               appears not to be current. Call lsetup with convfail=CV_FAIL_BAD_J\r
+               and then loop again */\r
+            convfail = CV_FAIL_BAD_J;\r
+\r
+            /* Rename some vectors for readibility */\r
+            vtemp1 = cv_mem.cv_tempv;\r
+            vtemp2 = cv_mem.cv_yS[0];\r
+            vtemp3 = cv_mem.cv_ftempS[0];\r
+\r
+            /* Call linear solver setup at converged y */\r
+            ier = cv_lsetup(cv_mem, convfail, cv_mem.cv_y,\r
+                                    cv_mem.cv_ftemp, cv_mem.cv_jcur, \r
+                                    vtemp1, vtemp2, vtemp3, cv_mem.cv_lsetup_select);\r
+            cv_mem.cv_nsetups++;\r
+            cv_mem.cv_nsetupsS++;\r
+            cv_mem.cv_gamrat = ONE;\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_gammap = cv_mem.cv_gamma;\r
+            cv_mem.cv_nstlp = cv_mem.cv_nst;\r
+\r
+            /* Return if lsetup failed */\r
+            if (ier < 0) return(CV_LSETUP_FAIL);\r
+            if (ier > 0) return(CONV_FAIL);\r
+\r
+          } /* end loop */\r
+\r
+        }\r
+\r
+        /*\r
+         * cvStgr1NewtonIteration\r
+         *\r
+         * This routine performs the Newton iteration for the i-th sensitivity. \r
+         * If the iteration succeeds, it returns the value CV_SUCCESS. \r
+         * If not, it may signal the cvStgr1NlsNewton routine to call lsetup \r
+         * and reattempt the iteration, by returning the value TRY_AGAIN. \r
+         * (In this case, cvStgr1NlsNewton must set convfail to CV_FAIL_BAD_J \r
+         * before calling setup again). Otherwise, this routine returns one \r
+         * of the appropriate values CV_LSOLVE_FAIL or CONV_FAIL back to \r
+         * cvStgr1NlsNewton.\r
+         */\r
+\r
+        private int cvStgr1NewtonIteration(CVodeMemRec cv_mem, int is)\r
+        {\r
+          int m, retval;\r
+          double Del, Delp, dcon;\r
+          NVector bS[];\r
+          NVector wrk1, wrk2;\r
+\r
+          m = 0;\r
+\r
+          /* Initialize Delp to avoid compiler warning message */\r
+          Delp = ZERO;\r
+\r
+          /* ftemp      <- f(t_n, y_n)\r
+             y          <- y_n\r
+             ftempS[is] <- fS(is, t_n, y_n(0), s_n(0))\r
+             acorS[is]  <- 0\r
+             yS[is]     <- yS_n(0)[is]                 */\r
+\r
+          for(;;) {\r
+\r
+            /* Evaluate the residual of the nonlinear systems */\r
+            N_VLinearSum_Serial(cv_mem.cv_rl1, cv_mem.cv_znS[1][is], ONE,\r
+                         cv_mem.cv_acorS[is], cv_mem.cv_tempvS[is]);\r
+            N_VLinearSum_Serial(cv_mem.cv_gamma, cv_mem.cv_ftempS[is], -ONE,\r
+                         cv_mem.cv_tempvS[is], cv_mem.cv_tempvS[is]);\r
+\r
+            /* Call the lsolve function */\r
+            bS = cv_mem.cv_tempvS;\r
+\r
+            cv_mem.cv_nniS1[is]++;\r
+\r
+            retval = cv_lsolve(cv_mem, bS[is], cv_mem.cv_ewtS[is],\r
+                                       cv_mem.cv_y, cv_mem.cv_ftemp,cv_mem.cv_lsolve_select);\r
+\r
+            /* Unrecoverable error in lsolve */\r
+            if (retval < 0) return(CV_LSOLVE_FAIL);\r
+\r
+            /* Recoverable error in lsolve and Jacobian data not current */\r
+            if (retval > 0) { \r
+              if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                return(TRY_AGAIN);\r
+              else\r
+                return(CONV_FAIL);\r
+            }\r
+\r
+            /* Get norm of correction; add correction to acorS and yS */\r
+            Del = N_VWrmsNorm_Serial(bS[is], cv_mem.cv_ewtS[is]);\r
+            N_VLinearSum_Serial(ONE, cv_mem.cv_acorS[is], ONE,\r
+                         bS[is], cv_mem.cv_acorS[is]);\r
+            N_VLinearSum_Serial(ONE, cv_mem.cv_znS[0][is], ONE,\r
+                         cv_mem.cv_acorS[is], cv_mem.cv_yS[is]);\r
+\r
+            /* Test for convergence.  If m > 0, an estimate of the convergence\r
+               rate constant is stored in crateS, and used in the test.        */\r
+            if (m > 0)\r
+              cv_mem.cv_crateS = Math.max(CRDOWN * cv_mem.cv_crateS, Del/Delp);\r
+            dcon = Del * Math.min(ONE, cv_mem.cv_crateS) / cv_mem.cv_tq[4];\r
+            if (dcon <= ONE) {\r
+              cv_mem.cv_jcur[0] = false;\r
+              return(CV_SUCCESS);  /* Convergence achieved */\r
+            }\r
+\r
+            m++;\r
+\r
+            /* Stop at maxcor iterations or if iter. seems to be diverging.\r
+               If still not converged and Jacobian data is not current, \r
+               signal to try the solution again                            */\r
+            if ((m == cv_mem.cv_maxcorS) || ((m >= 2) && (Del > RDIV * Delp))) {\r
+              if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                return(TRY_AGAIN);\r
+              else\r
+                return(CONV_FAIL);\r
+            }\r
+\r
+            /* Save norm of correction, evaluate fS, and loop again */\r
+            Delp = Del;\r
+\r
+            wrk1 = cv_mem.cv_tempv;\r
+            wrk2 = cv_mem.cv_tempvS[0];\r
+            retval = cvSensRhs1Wrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_y,\r
+                                       cv_mem.cv_ftemp, is, cv_mem.cv_yS[is], \r
+                                       cv_mem.cv_ftempS[is], wrk1, wrk2);\r
+            if (retval < 0) return(CV_SRHSFUNC_FAIL);\r
+            if (retval > 0) {\r
+              if ((!cv_mem.cv_jcur[0]) && (cv_mem.cv_lsetup_select > 0))\r
+                return(TRY_AGAIN);\r
+              else\r
+                return(SRHSFUNC_RECVR);\r
+            }\r
+\r
+          } /* end loop */\r
+\r
+        }\r
+        \r
+        /*\r
+         * cvSensRhs1Wrapper\r
+         *\r
+         * cvSensRhs1Wrapper is a high level routine that returns right-hand\r
+         * side of the is-th sensitivity equation. \r
+         *\r
+         * cvSensRhs1Wrapper is called only during the CV_STAGGERED1 corrector loop\r
+         * (ifS must be CV_ONESENS, otherwise CVodeSensInit would have \r
+         * issued an error message).\r
+         *\r
+         * The return value is that of the sensitivity RHS function fS1,\r
+         */\r
+\r
+        private int cvSensRhs1Wrapper(CVodeMemRec cv_mem, double time, \r
+                              NVector ycur, NVector fcur, \r
+                              int is, NVector yScur, NVector fScur,\r
+                              NVector temp1, NVector temp2)\r
+        {\r
+          int retval;\r
+\r
+          retval = cvSensRhs1InternalDQ(cv_mem.cv_Ns, time, ycur, fcur, is, yScur, \r
+                                  fScur, cv_mem.cv_fS_data, temp1, temp2);\r
+          cv_mem.cv_nfSe++;\r
+\r
+          return(retval);\r
+        }\r
+\r
+        /*\r
+         * cvSensRhsInternalDQ   - internal CVSensRhsFn\r
+         *\r
+         * cvSensRhsInternalDQ computes right hand side of all sensitivity equations\r
+         * by finite differences\r
+         */\r
+\r
+        private int cvSensRhsInternalDQ(int Ns, double t, \r
+                                NVector y, NVector ydot, \r
+                                NVector yS[], NVector ySdot[], \r
+                                CVodeMemRec cv_mem,  \r
+                                NVector ytemp, NVector ftemp)\r
+        {\r
+          int is, retval;\r
+          \r
+          for (is=0; is<Ns; is++) {\r
+            retval = cvSensRhs1InternalDQ(Ns, t, y, ydot, is, yS[is], \r
+                                          ySdot[is], cv_mem, ytemp, ftemp);\r
+            if (retval!=0) return(retval);\r
+          }\r
+\r
+          return(0);\r
+        }\r
+\r
+        /*\r
+         * cvSensRhs1InternalDQ   - internal CVSensRhs1Fn\r
+         *\r
+         * cvSensRhs1InternalDQ computes the right hand side of the is-th sensitivity \r
+         * equation by finite differences\r
+         *\r
+         * cvSensRhs1InternalDQ returns 0 if successful. Otherwise it returns the \r
+         * non-zero return value from f().\r
+         */\r
+\r
+        private int cvSensRhs1InternalDQ(int Ns, double t, \r
+                                 NVector y, NVector ydot, \r
+                                 int is, NVector yS, NVector ySdot, \r
+                                 CVodeMemRec cv_mem,\r
+                                 NVector ytemp, NVector ftemp)\r
+        {\r
+          int retval, method;\r
+          int nfel = 0, which;\r
+          double psave, pbari;\r
+          double delta , rdelta;\r
+          double Deltap, rDeltap, r2Deltap;\r
+          double Deltay, rDeltay, r2Deltay;\r
+          double Delta , rDelta , r2Delta ;\r
+          double norms, ratio;\r
+          \r
+          /* cvode_mem is passed here as user data */\r
+\r
+          delta = Math.sqrt(Math.max(cv_mem.cv_reltol, cv_mem.cv_uround));\r
+          rdelta = ONE/delta;\r
+          \r
+          pbari = cv_mem.cv_pbar[is];\r
+\r
+          which = cv_mem.cv_plist[is];\r
+\r
+          psave = cv_mem.cv_p[which];\r
+          \r
+          Deltap  = pbari * delta;\r
+          rDeltap = ONE/Deltap;\r
+          norms   = N_VWrmsNorm_Serial(yS, cv_mem.cv_ewt) * pbari;\r
+          rDeltay = Math.max(norms, rdelta) / pbari;\r
+          Deltay  = ONE/rDeltay;\r
+          \r
+          if (cv_mem.cv_DQrhomax == ZERO) {\r
+            /* No switching */\r
+            method = (cv_mem.cv_DQtype==CV_CENTERED) ? CENTERED1 : FORWARD1;\r
+          } else {\r
+            /* switch between simultaneous/separate DQ */\r
+            ratio = Deltay * rDeltap;\r
+            if ( Math.max(ONE/ratio, ratio) <= cv_mem.cv_DQrhomax )\r
+              method = (cv_mem.cv_DQtype==CV_CENTERED) ? CENTERED1 : FORWARD1;\r
+            else\r
+              method = (cv_mem.cv_DQtype==CV_CENTERED) ? CENTERED2 : FORWARD2;\r
+          }\r
+\r
+          switch(method) {\r
+            \r
+          case CENTERED1:\r
+            \r
+            Delta = Math.min(Deltay, Deltap);\r
+            r2Delta = HALF/Delta;\r
+            \r
+            N_VLinearSum_Serial(ONE,y,Delta,yS,ytemp);\r
+            cv_mem.cv_p[which] = psave + Delta;\r
+\r
+            retval = f(t, ytemp, ySdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+            \r
+            N_VLinearSum_Serial(ONE,y,-Delta,yS,ytemp);\r
+            cv_mem.cv_p[which] = psave - Delta;\r
+\r
+            retval = f(t, ytemp, ftemp, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(r2Delta,ySdot,-r2Delta,ftemp,ySdot);\r
+            \r
+            break;\r
+            \r
+          case CENTERED2:\r
+            \r
+            r2Deltap = HALF/Deltap;\r
+            r2Deltay = HALF/Deltay;\r
+            \r
+            N_VLinearSum_Serial(ONE,y,Deltay,yS,ytemp);\r
+\r
+            retval = f(t, ytemp, ySdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(ONE,y,-Deltay,yS,ytemp);\r
+\r
+            retval = f(t, ytemp, ftemp, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(r2Deltay, ySdot, -r2Deltay, ftemp, ySdot);\r
+            \r
+            cv_mem.cv_p[which] = psave + Deltap;\r
+            retval = f(t, y, ytemp, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            cv_mem.cv_p[which] = psave - Deltap;\r
+            retval = f(t, y, ftemp, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(r2Deltap,ytemp,-r2Deltap,ftemp,ftemp);\r
+            \r
+            N_VLinearSum_Serial(ONE,ySdot,ONE,ftemp,ySdot);\r
+            \r
+            break;\r
+            \r
+          case FORWARD1:\r
+            \r
+            Delta = Math.min(Deltay, Deltap);\r
+            rDelta = ONE/Delta;\r
+            \r
+            N_VLinearSum_Serial(ONE,y,Delta,yS,ytemp);\r
+            cv_mem.cv_p[which] = psave + Delta;\r
+\r
+            retval = f(t, ytemp, ySdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+            \r
+            N_VLinearSum_Serial(rDelta,ySdot,-rDelta,ydot,ySdot);\r
+            \r
+            break;\r
+            \r
+          case FORWARD2:\r
+            \r
+            N_VLinearSum_Serial(ONE,y,Deltay,yS,ytemp);\r
+\r
+            retval = f(t, ytemp, ySdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(rDeltay, ySdot, -rDeltay, ydot, ySdot);\r
+            \r
+            cv_mem.cv_p[which] = psave + Deltap;\r
+            retval = f(t, y, ytemp, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(rDeltap,ytemp,-rDeltap,ydot,ftemp);\r
+              \r
+            N_VLinearSum_Serial(ONE,ySdot,ONE,ftemp,ySdot);\r
+            \r
+            break;\r
+            \r
+          }\r
+          \r
+          cv_mem.cv_p[which] = psave;\r
+          \r
+          /* Increment counter nfeS */\r
+          cv_mem.cv_nfeS += nfel;\r
+          \r
+          return(0);\r
+        }\r
+        \r
+        /*\r
+         * cvQuadSensRhsInternalDQ   - internal CVQuadSensRhsFn\r
+         *\r
+         * cvQuadSensRhsInternalDQ computes right hand side of all quadrature\r
+         * sensitivity equations by finite differences. All work is actually\r
+         * done in cvQuadSensRhs1InternalDQ.\r
+         */\r
+\r
+        private int cvQuadSensRhsInternalDQ(int Ns, double t, \r
+                                           NVector y, NVector yS[],\r
+                                           NVector yQdot, NVector yQSdot[],\r
+                                           CVodeMemRec cv_mem,  \r
+                                           NVector tmp, NVector tmpQ)\r
+        {\r
+          int is, retval;\r
+          \r
+          /* cvode_mem is passed here as user data */\r
+\r
+          for (is=0; is<Ns; is++) {\r
+            retval = cvQuadSensRhs1InternalDQ(cv_mem, is, t,\r
+                                              y, yS[is], \r
+                                              yQdot, yQSdot[is],\r
+                                              tmp, tmpQ);\r
+            if (retval!=0) return(retval);\r
+          }\r
+\r
+          return(0);\r
+        }\r
+\r
+        private int cvQuadSensRhs1InternalDQ(CVodeMemRec cv_mem, int is, double t, \r
+                                            NVector y, NVector yS,\r
+                                            NVector yQdot, NVector yQSdot, \r
+                                            NVector tmp, NVector tmpQ)\r
+        {\r
+          int retval, method;\r
+          int nfel = 0, which;\r
+          double psave, pbari;\r
+          double delta , rdelta;\r
+          double Deltap;\r
+          double Deltay, rDeltay;\r
+          double Delta , rDelta , r2Delta ;\r
+          double norms;\r
+\r
+          delta = Math.sqrt(Math.max(cv_mem.cv_reltol, cv_mem.cv_uround));\r
+          rdelta = ONE/delta;\r
+          \r
+          pbari = cv_mem.cv_pbar[is];\r
+\r
+          which = cv_mem.cv_plist[is];\r
+\r
+          psave = cv_mem.cv_p[which];\r
+          \r
+          Deltap  = pbari * delta;\r
+          norms   = N_VWrmsNorm_Serial(yS, cv_mem.cv_ewt) * pbari;\r
+          rDeltay = Math.max(norms, rdelta) / pbari;\r
+          Deltay  = ONE/rDeltay;\r
+          \r
+          method = (cv_mem.cv_DQtype==CV_CENTERED) ? CENTERED1 : FORWARD1;\r
+\r
+          switch(method) {\r
+\r
+          case CENTERED1:\r
+            \r
+            Delta = Math.min(Deltay, Deltap);\r
+            r2Delta = HALF/Delta;\r
+            \r
+            N_VLinearSum_Serial(ONE, y, Delta, yS, tmp);\r
+            cv_mem.cv_p[which] = psave + Delta;\r
+\r
+            retval = fQ(t, tmp, yQSdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+            \r
+            N_VLinearSum_Serial(ONE, y, -Delta, yS, tmp);\r
+            cv_mem.cv_p[which] = psave - Delta;\r
+\r
+            retval = fQ(t, tmp, tmpQ, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+\r
+            N_VLinearSum_Serial(r2Delta, yQSdot, -r2Delta, tmpQ, yQSdot);\r
+            \r
+            break;\r
+\r
+          case FORWARD1:\r
+            \r
+            Delta = Math.min(Deltay, Deltap);\r
+            rDelta = ONE/Delta;\r
+            \r
+            N_VLinearSum_Serial(ONE, y, Delta, yS, tmp);\r
+            cv_mem.cv_p[which] = psave + Delta;\r
+\r
+            retval = fQ(t, tmp, yQSdot, cv_mem.cv_user_data);\r
+            nfel++;\r
+            if (retval != 0) return(retval);\r
+            \r
+            N_VLinearSum_Serial(rDelta, yQSdot, -rDelta, yQdot, yQSdot);\r
+            \r
+            break;\r
+\r
+          }\r
+\r
+          cv_mem.cv_p[which] = psave;\r
+          \r
+          /* Increment counter nfQeS */\r
+          cv_mem.cv_nfQeS += nfel;\r
+          \r
+          return(0);\r
+        }\r
 \r
 \r
 }
\ No newline at end of file