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

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

index fbbc79f..1b809b2 100644 (file)
@@ -593,10 +593,14 @@ public abstract class CVODES {
     // For cv_mem.cv_setup_select\r
     final int cvDlsSetup_select = 1;\r
     final int cvDlsLapackSetup_select = 2;\r
-    // For cv_mem.cv_fS1_select\r
+    // For cv_mem.cv_fS1\r
     final int cvSensRhs1InternalDQ_select = -1;\r
-    // For cv_mem.cv_f_select\r
+    // For cv_mem.cv_f\r
     final int cvSensRhsInternalDQ_select = -1;\r
+    final int cvDlsJacBWrapper_select = 101;\r
+    final int CVArhs_select = 100;\r
+       final int cvDlsFreeB_select = 100;\r
+       final int CVArhsQ_select = 200;\r
 \r
     final int cvsRoberts_dns = 1;\r
     final int cvsDirectDemo_ls_Problem_1 = 2;\r
@@ -2530,9 +2534,9 @@ public abstract class CVODES {
          int cv_ism;                 /* ism = SIMULTANEOUS or STAGGERED              */\r
 \r
          //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
-         int cv_fS_select;\r
+         int cv_fS;\r
          //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
-         int cv_fS1_select;\r
+         int cv_fS1;\r
          //void *cv_fS_data;           /* data pointer passed to fS                    */\r
          UserData cv_fS_data = new UserData();;\r
          boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
@@ -2966,7 +2970,7 @@ public abstract class CVODES {
          /* Set default values for quad. optional inputs */\r
 \r
          cv_mem.cv_quadr      = false;\r
-         //cv_mem.cv_fQ         = null;\r
+         cv_mem.cv_fQ         = -1;\r
          cv_mem.cv_errconQ    = false;\r
          cv_mem.cv_itolQ      = CV_NN;\r
 \r
@@ -2975,9 +2979,9 @@ public abstract class CVODES {
          cv_mem.cv_sensi      = false;\r
          cv_mem.cv_fS_data    = new UserData();\r
          //cv_mem.cv_fS         = cvSensRhsInternalDQ;\r
-         cv_mem.cv_fS_select = cvSensRhsInternalDQ_select;\r
+         cv_mem.cv_fS = cvSensRhsInternalDQ_select;\r
          //cv_mem.cv_fS1         = cvSensRhs1InternalDQ;\r
-         cv_mem.cv_fS1_select = cvSensRhs1InternalDQ_select;\r
+         cv_mem.cv_fS1 = cvSensRhs1InternalDQ_select;\r
          cv_mem.cv_fSDQ       = true;\r
          cv_mem.cv_ifS        = CV_ONESENS;\r
          cv_mem.cv_DQtype     = CV_CENTERED;\r
@@ -3139,7 +3143,7 @@ public abstract class CVODES {
 \r
        /* Copy the input parameters into CVODES state */\r
 \r
-       //cv_mem.cv_f  = f;\r
+       cv_mem.cv_f  = f;\r
        cv_mem.cv_tn = t0;\r
 \r
        /* Set step parameters */\r
@@ -3804,6 +3808,7 @@ public abstract class CVODES {
     }\r
     cvdls_mem = cv_mem.cv_lmem;\r
 \r
+    \r
     if (jac >= 0) {\r
       cvdls_mem.jacDQ  = false;\r
       cvdls_mem.jac    = jac;\r
@@ -3914,7 +3919,11 @@ public abstract class CVODES {
        * If computing quadr. sensi., call fQS at (t0,y0,yS0), set znQS[1][is] = yQS'(t0), is=1,...,Ns.\r
        */\r
 \r
-      if (testMode) {\r
+      if (cv_mem.cv_f == CVArhs_select) {\r
+         retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                  cv_mem.cv_zn[1], cv_mem.cv_user_data.memRec);     \r
+      }\r
+      else if (testMode) {\r
          retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                   cv_mem.cv_zn[1], cv_mem.cv_user_data);   \r
       }\r
@@ -3935,7 +3944,11 @@ public abstract class CVODES {
       }\r
 \r
       if (cv_mem.cv_quadr) {\r
-       if (testMode) {\r
+       if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+               retval = CVArhsQ(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                    cv_mem.cv_znQ[1], cv_mem.cv_user_data.memRec);             \r
+       }\r
+        else if (testMode) {\r
                retval = fQTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                     cv_mem.cv_znQ[1], cv_mem.cv_user_data);    \r
        }\r
@@ -5801,7 +5814,11 @@ public abstract class CVODES {
    \r
    /* tempv <- f(t+h, h*y'(t)+y(t)) */\r
 \r
-   if (testMode) {\r
+   if (cv_mem.cv_f == CVArhs_select) {\r
+          retval = CVArhs(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
+               cv_mem.cv_tempv, cv_mem.cv_user_data.memRec);      \r
+   }\r
+   else if (testMode) {\r
           retval = fTestMode(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
                cv_mem.cv_tempv, cv_mem.cv_user_data);   \r
    }\r
@@ -5814,7 +5831,11 @@ public abstract class CVODES {
    if (retval > 0) return(RHSFUNC_RECVR);\r
 \r
    if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
-        if (testMode) {\r
+        if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+                retval = CVArhsQ(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
+                 cv_mem.cv_tempvQ, cv_mem.cv_user_data.memRec);         \r
+        }\r
+     else if (testMode) {\r
                 retval = fQTestMode(cv_mem.cv_tn+hg, cv_mem.cv_y,\r
                  cv_mem.cv_tempvQ, cv_mem.cv_user_data);        \r
         }\r
@@ -6754,8 +6775,11 @@ else                return(snrm);
        /* 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
-       if (testMode) {\r
+       if (cv_mem.cv_f == CVArhs_select) {\r
+          retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_y,\r
+                   cv_mem.cv_ftemp, cv_mem.cv_user_data.memRec);      \r
+       }\r
+       else if (testMode) {\r
           retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_y,\r
                    cv_mem.cv_ftemp, cv_mem.cv_user_data);   \r
        }\r
@@ -7558,7 +7582,11 @@ else                return(snrm);
    /* Initialize delS and Delp to avoid compiler warning message */\r
    delS = Delp = ZERO;\r
 \r
-   if (testMode) {\r
+   if (cv_mem.cv_f == CVArhs_select) {\r
+          retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+               cv_mem.cv_tempv, cv_mem.cv_user_data.memRec);    \r
+   }\r
+   else if (testMode) {\r
           retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                cv_mem.cv_tempv, cv_mem.cv_user_data);   \r
    }\r
@@ -7666,7 +7694,11 @@ else                return(snrm);
 \r
      Delp = Del;\r
 \r
-     if (testMode) {\r
+     if (cv_mem.cv_f == CVArhs_select) {\r
+        retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_y,\r
+                 cv_mem.cv_tempv, cv_mem.cv_user_data.memRec);  \r
+     }\r
+     else if (testMode) {\r
         retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_y,\r
                  cv_mem.cv_tempv, cv_mem.cv_user_data);\r
         \r
@@ -7768,7 +7800,11 @@ else                return(snrm);
 \r
    for(;;) {\r
 \r
-     if (testMode) {\r
+     if (cv_mem.cv_f == CVArhs_select) {\r
+        retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                 cv_mem.cv_ftemp, cv_mem.cv_user_data.memRec);          \r
+     }\r
+        else if (testMode) {\r
         retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                  cv_mem.cv_ftemp, cv_mem.cv_user_data);         \r
      }\r
@@ -7912,6 +7948,10 @@ else                return(snrm);
        if (cvdls_mem.jacDQ) {\r
           cvDlsDenseDQJac(cv_mem.cv_tn, y, fy, cvdls_mem.A, cv_mem, tmp1);  \r
        }\r
+       else if (cvdls_mem.jac == cvDlsJacBWrapper_select) {\r
+          retval = cvDlsJacBWrapper(cv_mem.cv_tn, y, fy, cvdls_mem.A, \r
+                   cvdls_mem.J_data.memRec, tmp1, tmp2, tmp3);  \r
+       }\r
        else if (testMode) {\r
           retval = JacTestMode(cv_mem.cv_tn, y, fy, cvdls_mem.A, \r
                    cvdls_mem.J_data, tmp1, tmp2, tmp3);   \r
@@ -8293,7 +8333,11 @@ else                return(snrm);
        \r
        /* Save norm of correction, evaluate f, and loop again */\r
        Delp = Del;\r
-       if (testMode) {\r
+       if (cv_mem.cv_f == CVArhs_select) {\r
+          retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_y,\r
+                   cv_mem.cv_ftemp, cv_mem.cv_user_data.memRec);      \r
+       }\r
+       else if (testMode) {\r
           retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_y,\r
                    cv_mem.cv_ftemp, cv_mem.cv_user_data);   \r
        }\r
@@ -8732,7 +8776,11 @@ else                return(snrm);
           cv_mem.cv_qwait = LONG_WAIT;\r
           cv_mem.cv_nscon = 0;\r
 \r
-          if (testMode) {\r
+          if (cv_mem.cv_f == CVArhs_select) {\r
+                  retval = CVArhs(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                   cv_mem.cv_tempv, cv_mem.cv_user_data.memRec);         \r
+          }\r
+          else if (testMode) {\r
                   retval = fTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                    cv_mem.cv_tempv, cv_mem.cv_user_data);   \r
           }\r
@@ -8748,7 +8796,11 @@ else                return(snrm);
 \r
           if (cv_mem.cv_quadr) {\r
 \r
-            if (testMode) {\r
+                if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+                        retval = CVArhsQ(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
+                     cv_mem.cv_tempvQ, cv_mem.cv_user_data.memRec);               \r
+                }\r
+            else if (testMode) {\r
                 retval = fQTestMode(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
                      cv_mem.cv_tempvQ, cv_mem.cv_user_data);    \r
             }\r
@@ -8823,7 +8875,11 @@ else                return(snrm);
           int retval;\r
 \r
           /* Save quadrature correction in acorQ */\r
-          if (testMode) {\r
+       if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+          retval = CVArhsQ(cv_mem.cv_tn, cv_mem.cv_y,\r
+                   cv_mem.cv_acorQ, cv_mem.cv_user_data.memRec);   \r
+       }\r
+          else if (testMode) {\r
                   retval = fQTestMode(cv_mem.cv_tn, cv_mem.cv_y,\r
                    cv_mem.cv_acorQ, cv_mem.cv_user_data);   \r
           }\r
@@ -9598,7 +9654,10 @@ else                return(snrm);
             N_VLinearSum_Serial(ONE,y,Delta,yS,ytemp);\r
             cv_mem.cv_p[which] = psave + Delta;\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ySdot, cv_mem.cv_user_data.memRec);           \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ySdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9610,7 +9669,10 @@ else                return(snrm);
             N_VLinearSum_Serial(ONE,y,-Delta,yS,ytemp);\r
             cv_mem.cv_p[which] = psave - Delta;\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ftemp, cv_mem.cv_user_data.memRec);   \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ftemp, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9630,7 +9692,10 @@ else                return(snrm);
             \r
             N_VLinearSum_Serial(ONE,y,Deltay,yS,ytemp);\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ySdot, cv_mem.cv_user_data.memRec);           \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ySdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9641,7 +9706,10 @@ else                return(snrm);
 \r
             N_VLinearSum_Serial(ONE,y,-Deltay,yS,ytemp);\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ftemp, cv_mem.cv_user_data.memRec);   \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ftemp, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9653,7 +9721,10 @@ else                return(snrm);
             N_VLinearSum_Serial(r2Deltay, ySdot, -r2Deltay, ftemp, ySdot);\r
             \r
             cv_mem.cv_p[which] = psave + Deltap;\r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, y, ytemp, cv_mem.cv_user_data.memRec);       \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, y, ytemp, cv_mem.cv_user_data);   \r
             }\r
             else {\r
@@ -9663,7 +9734,10 @@ else                return(snrm);
             if (retval != 0) return(retval);\r
 \r
             cv_mem.cv_p[which] = psave - Deltap;\r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, y, ftemp, cv_mem.cv_user_data.memRec);               \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, y, ftemp, cv_mem.cv_user_data);   \r
             }\r
             else {\r
@@ -9686,7 +9760,10 @@ else                return(snrm);
             N_VLinearSum_Serial(ONE,y,Delta,yS,ytemp);\r
             cv_mem.cv_p[which] = psave + Delta;\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ySdot, cv_mem.cv_user_data.memRec);   \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ySdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9703,7 +9780,10 @@ else                return(snrm);
             \r
             N_VLinearSum_Serial(ONE,y,Deltay,yS,ytemp);\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, ytemp, ySdot, cv_mem.cv_user_data.memRec);           \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, ytemp, ySdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9715,7 +9795,10 @@ else                return(snrm);
             N_VLinearSum_Serial(rDeltay, ySdot, -rDeltay, ydot, ySdot);\r
             \r
             cv_mem.cv_p[which] = psave + Deltap;\r
-            if (testMode) {\r
+            if (cv_mem.cv_f == CVArhs_select) {\r
+                retval = CVArhs(t, y, ytemp, cv_mem.cv_user_data.memRec);               \r
+            }\r
+            else if (testMode) {\r
                 retval = fTestMode(t, y, ytemp, cv_mem.cv_user_data);   \r
             }\r
             else {\r
@@ -9809,7 +9892,10 @@ else                return(snrm);
             N_VLinearSum_Serial(ONE, y, Delta, yS, tmp);\r
             cv_mem.cv_p[which] = psave + Delta;\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+                retval = CVArhsQ(t, tmp, yQSdot, cv_mem.cv_user_data.memRec);           \r
+            }\r
+            else if (testMode) {\r
                 retval = fQTestMode(t, tmp, yQSdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -9820,8 +9906,11 @@ else                return(snrm);
             \r
             N_VLinearSum_Serial(ONE, y, -Delta, yS, tmp);\r
             cv_mem.cv_p[which] = psave - Delta;\r
-\r
-            if (testMode) {\r
+         \r
+            if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+                retval = CVArhsQ(t, tmp, tmpQ, cv_mem.cv_user_data.memRec);     \r
+            }\r
+            else if (testMode) {\r
                 retval = fQTestMode(t, tmp, tmpQ, cv_mem.cv_user_data);         \r
             }\r
             else {\r
@@ -9842,7 +9931,10 @@ else                return(snrm);
             N_VLinearSum_Serial(ONE, y, Delta, yS, tmp);\r
             cv_mem.cv_p[which] = psave + Delta;\r
 \r
-            if (testMode) {\r
+            if (cv_mem.cv_fQ == CVArhsQ_select) {\r
+                retval = CVArhsQ(t, tmp, yQSdot, cv_mem.cv_user_data.memRec);           \r
+            }\r
+            else if (testMode) {\r
                 retval = fQTestMode(t, tmp, yQSdot, cv_mem.cv_user_data);       \r
             }\r
             else {\r
@@ -10948,7 +11040,7 @@ else                return(snrm);
             return(CV_SUCCESS);\r
           }\r
 \r
-          private int CVodeSStolerances(CVodeMemRec cv_mem, double reltol, double abstol)\r
+          protected int CVodeSStolerances(CVodeMemRec cv_mem, double reltol, double abstol)\r
           {\r
 \r
             if (cv_mem==null) {\r
@@ -11392,18 +11484,18 @@ else                return(snrm);
 \r
          cv_mem.cv_ifS = CV_ONESENS;\r
          //cv_mem.cv_fS  = null;\r
-         cv_mem.cv_fS_select = -1;\r
+         cv_mem.cv_fS = -1;\r
 \r
          if (fS1_select < 0) {\r
 \r
            cv_mem.cv_fSDQ = true;\r
-           cv_mem.cv_fS1_select  = cvSensRhs1InternalDQ_select;\r
+           cv_mem.cv_fS1  = cvSensRhs1InternalDQ_select;\r
            cv_mem.cv_fS_data.memRec = cv_mem;\r
 \r
          } else {\r
 \r
            cv_mem.cv_fSDQ = false;\r
-           cv_mem.cv_fS1_select  = fS1_select;\r
+           cv_mem.cv_fS1  = fS1_select;\r
            cv_mem.cv_fS_data = cv_mem.cv_user_data;\r
 \r
          }\r
@@ -11924,20 +12016,23 @@ else                return(snrm);
 \r
                  /* Right hand side function for backward run */\r
                  //CVRhsFnB cv_f;\r
+                 int cv_f;\r
                  //CVRhsFnBS cv_fs;\r
 \r
                  /* Right hand side quadrature function for backward run */\r
                  //CVQuadRhsFnB cv_fQ;\r
+                 int cv_fQ;\r
                  //CVQuadRhsFnBS cv_fQs;\r
 \r
                  /* User user_data */\r
                  UserData cv_user_data;\r
                    \r
                  /* Memory block for a linear solver's interface to CVODEA */\r
-                 //void *cv_lmem;\r
+                 CVDlsMemRecB cv_lmem;\r
 \r
                  /* Function to free any memory allocated by the linear solver */\r
                  //int (*cv_lfree)(CVodeBMem cvB_mem);\r
+                 int cv_lfree_select;\r
 \r
                  /* Memory block for a preconditioner's module interface to CVODEA */ \r
                  //void *cv_pmem;\r
@@ -12058,5 +12153,509 @@ else                return(snrm);
                  int order;\r
                };\r
 \r
+               class CVDlsMemRecB {\r
+\r
+                         //CVDlsJacFnB jacB;\r
+                         int jacB;\r
+                         //CVDlsJacFnBS jacBS;\r
+                 int jacBS;\r
+                       };\r
+                       \r
+               /* cvDlsJacBWrapper interfaces to the CVDlsJacFnB routine provided \r
+                  by the user. cvDlsJacBWrapper is of type CVDlsJacFn.\r
+                  NOTE: data here contains cvode_mem */\r
+               private int cvDlsJacBWrapper(double t, NVector yB, NVector fyB, \r
+                                           double JB[][], CVodeMemRec cv_mem,\r
+                                           NVector tmp1B, NVector tmp2B, NVector tmp3B)\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 CVodeBMemRec cvB_mem;\r
+                 int retval, flag;\r
+\r
+                 /* Check if cvode_mem exists */\r
+                 if (cv_mem == null) {\r
+                   cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS",\r
+                                  "cvDlsJacBWrapper", MSGD_CVMEM_NULL);\r
+                   return(CVDLS_MEM_NULL);\r
+                 }\r
+\r
+                 /* Was ASA initialized? */\r
+                 if (cv_mem.cv_adjMallocDone == false) {\r
+                   cvProcessError(cv_mem, CVDLS_NO_ADJ, "CVSDLS",\r
+                                  "cvDlsJacBWrapper", MSGD_NO_ADJ);\r
+                   return(CVDLS_NO_ADJ);\r
+                 } \r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                 if (ca_mem.ca_bckpbCrt == null) {\r
+                   cvProcessError(cv_mem, CVDLS_LMEMB_NULL, "CVSDLS",\r
+                                  "cvDlsJacBWrapper", MSGD_LMEMB_NULL);\r
+                   return(CVDLS_LMEMB_NULL);\r
+                 }\r
+                 cvB_mem = ca_mem.ca_bckpbCrt;\r
+\r
+                 /* Forward solution from interpolation */\r
+                 if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                       flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                 }\r
+                 else {\r
+                         flag = CVApolynomialGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                 }\r
+                 if (flag != CV_SUCCESS) {\r
+                   cvProcessError(cv_mem, -1, "CVSDLS", "cvDlsJacBWrapper",\r
+                                  MSGD_BAD_TINTERP);\r
+                   return(-1);\r
+                 }\r
+\r
+                 /* Call user's adjoint jacB routine (of type CVDlsJacFnB) */\r
+                 if (testMode) {\r
+                         retval = JacBTestMode(t, ca_mem.ca_ytmp, yB, fyB,\r
+                      JB, cvB_mem.cv_user_data,\r
+                      tmp1B, tmp2B, tmp3B);\r
+                 }\r
+                 else {\r
+                 retval = JacB(t, ca_mem.ca_ytmp, yB, fyB,\r
+                                           JB, cvB_mem.cv_user_data,\r
+                                           tmp1B, tmp2B, tmp3B);\r
+                 }\r
+                 return(retval);\r
+               }\r
+               \r
+               /*\r
+                * CVAhermiteGetY ( -> IMget )\r
+                *\r
+                * This routine uses cubic piece-wise Hermite interpolation for \r
+                * the forward solution vector. \r
+                * It is typically called by the wrapper routines before calling\r
+                * user provided routines (fB, djacB, bjacB, jtimesB, psolB) but\r
+                * can be directly called by the user through CVodeGetAdjY\r
+                */\r
+\r
+               protected int CVAhermiteGetY(CVodeMemRec cv_mem, double t,\r
+                                         NVector y, NVector yS[])\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 DtpntMemRec dt_mem[];\r
+                 HermiteDataMemRec content0, content1;\r
+\r
+                 double t0, t1, delta;\r
+                 double factor1, factor2, factor3;\r
+\r
+                 NVector y0, yd0, y1, yd1;\r
+                 NVector yS0[]=null;\r
+                 NVector ySd0[]= null;\r
+                 NVector yS1[], ySd1[];\r
+\r
+                 int flag, is, NS;\r
+                 long indx[] = new long[1];\r
+                 boolean newpoint[] = new boolean[1];\r
+\r
+                \r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+                 dt_mem = ca_mem.dt_mem;\r
+                \r
+                 /* Local value of Ns */\r
+                \r
+                 NS = (ca_mem.ca_IMinterpSensi && (yS != null)) ? cv_mem.cv_Ns : 0;\r
+\r
+                 /* Get the index in dt_mem */\r
+\r
+                 flag = CVAfindIndex(cv_mem, t, indx, newpoint);\r
+                 if (flag != CV_SUCCESS) return(flag);\r
+\r
+                 /* If we are beyond the left limit but close enough,\r
+                    then return y at the left limit. */\r
+\r
+                 if (indx[0] == 0) {\r
+                   content0 = dt_mem[0].hermiteContent;\r
+                   N_VScale_Serial(ONE, content0.y, y);\r
+                   for (is=0; is<NS; is++) N_VScale_Serial(ONE, content0.yS[is], yS[is]);\r
+                   return(CV_SUCCESS);\r
+                 }\r
+\r
+                 /* Extract stuff from the appropriate data points */\r
+\r
+                 t0 = dt_mem[(int)(indx[0]-1)].t;\r
+                 t1 = dt_mem[(int)indx[0]].t;\r
+                 delta = t1 - t0;\r
+\r
+                 content0 = dt_mem[(int)(indx[0]-1)].hermiteContent;\r
+                 y0  = content0.y;\r
+                 yd0 = content0.yd;\r
+                 if (ca_mem.ca_IMinterpSensi) {\r
+                   yS0  = content0.yS;\r
+                   ySd0 = content0.ySd;\r
+                 }\r
+\r
+                 if (newpoint[0]) {\r
+                   \r
+                   /* Recompute Y0 and Y1 */\r
+\r
+                   content1 = dt_mem[(int)indx[0]].hermiteContent;\r
+\r
+                   y1  = content1.y;\r
+                   yd1 = content1.yd;\r
+\r
+                   N_VLinearSum_Serial(ONE, y1, -ONE, y0, ca_mem.ca_Y[0]);\r
+                   N_VLinearSum_Serial(ONE, yd1,  ONE, yd0, ca_mem.ca_Y[1]);\r
+                   N_VLinearSum_Serial(delta, ca_mem.ca_Y[1], -TWO, ca_mem.ca_Y[0], ca_mem.ca_Y[1]);\r
+                   N_VLinearSum_Serial(ONE, ca_mem.ca_Y[0], -delta, yd0, ca_mem.ca_Y[0]);\r
+\r
+\r
+                   yS1  = content1.yS;\r
+                   ySd1 = content1.ySd;\r
+                     \r
+                   for (is=0; is<NS; is++) {\r
+                     N_VLinearSum_Serial(ONE, yS1[is], -ONE, yS0[is], ca_mem.ca_YS[0][is]);\r
+                     N_VLinearSum_Serial(ONE, ySd1[is],  ONE, ySd0[is], ca_mem.ca_YS[1][is]);\r
+                     N_VLinearSum_Serial(delta, ca_mem.ca_YS[1][is], -TWO, ca_mem.ca_YS[0][is], ca_mem.ca_YS[1][is]);\r
+                     N_VLinearSum_Serial(ONE, ca_mem.ca_YS[0][is], -delta, ySd0[is], ca_mem.ca_YS[0][is]);\r
+                   }\r
+\r
+                 }\r
+\r
+                 /* Perform the actual interpolation. */\r
+\r
+                 factor1 = t - t0;\r
+\r
+                 factor2 = factor1/delta;\r
+                 factor2 = factor2*factor2;\r
+\r
+                 factor3 = factor2*(t-t1)/delta;\r
+\r
+                 N_VLinearSum_Serial(ONE, y0, factor1, yd0, y);\r
+                 N_VLinearSum_Serial(ONE, y, factor2, ca_mem.ca_Y[0], y);\r
+                 N_VLinearSum_Serial(ONE, y, factor3, ca_mem.ca_Y[1], y);\r
+\r
+                 for (is=0; is<NS; is++) {\r
+                   N_VLinearSum_Serial(ONE, yS0[is], factor1, ySd0[is], yS[is]);\r
+                   N_VLinearSum_Serial(ONE, yS[is], factor2, ca_mem.ca_YS[0][is], yS[is]);\r
+                   N_VLinearSum_Serial(ONE, yS[is], factor3, ca_mem.ca_YS[1][is], yS[is]);\r
+                 }\r
+\r
+\r
+                 return(CV_SUCCESS);\r
+               }\r
+               \r
+               /*\r
+                * CVAfindIndex\r
+                *\r
+                * Finds the index in the array of data point strctures such that\r
+                *     dt_mem[indx-1].t <= t < dt_mem[indx].t\r
+                * If indx is changed from the previous invocation, then newpoint = SUNTRUE\r
+                *\r
+                * If t is beyond the leftmost limit, but close enough, indx=0.\r
+                *\r
+                * Returns CV_SUCCESS if successful and CV_GETY_BADT if unable to\r
+                * find indx (t is too far beyond limits).\r
+                */\r
+\r
+               private int CVAfindIndex(CVodeMemRec cv_mem, double t, \r
+                                       long indx[], boolean newpoint[])\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 long ilast = 0;\r
+                 DtpntMemRec dt_mem[];\r
+                 int sign;\r
+                 boolean to_left, to_right;\r
+\r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+                 dt_mem = ca_mem.dt_mem;\r
+\r
+                 newpoint[0] = false;\r
+\r
+                 /* Find the direction of integration */\r
+                 sign = (ca_mem.ca_tfinal - ca_mem.ca_tinitial > ZERO) ? 1 : -1;\r
+\r
+                 /* If this is the first time we use new data */\r
+                 if (ca_mem.ca_IMnewData) {\r
+                   ilast     = ca_mem.ca_np-1;\r
+                   newpoint[0] = true;\r
+                   ca_mem.ca_IMnewData   = false;\r
+                 }\r
+\r
+                 /* Search for indx starting from ilast */\r
+                 to_left  = ( sign*(t - dt_mem[(int)(ilast-1)].t) < ZERO);\r
+                 to_right = ( sign*(t - dt_mem[(int)(ilast)].t)   > ZERO);\r
+\r
+                 if ( to_left ) {\r
+                   /* look for a new indx to the left */\r
+\r
+                   newpoint[0] = true;\r
+                   \r
+                   indx[0] = ilast;\r
+                   for(;;) {\r
+                     if ( indx[0] == 0 ) break;\r
+                     if ( sign*(t - dt_mem[(int)(indx[0]-1)].t) <= ZERO ) (indx[0])--;\r
+                     else                                         break;\r
+                   }\r
+\r
+                   if ( indx[0] == 0 )\r
+                     ilast = 1;\r
+                   else\r
+                     ilast = indx[0];\r
+\r
+                   if ( indx[0] == 0 ) {\r
+                     /* t is beyond leftmost limit. Is it too far? */  \r
+                     if ( Math.abs(t - dt_mem[0].t) > FUZZ_FACTOR * cv_mem.cv_uround ) {\r
+                       return(CV_GETY_BADT);\r
+                     }\r
+                   }\r
+\r
+                 } else if ( to_right ) {\r
+                   /* look for a new indx to the right */\r
+\r
+                   newpoint[0] = true;\r
+\r
+                   indx [0]= ilast;\r
+                   for(;;) {\r
+                     if ( sign*(t - dt_mem[(int)indx[0]].t) > ZERO) (indx[0])++;\r
+                     else                                     break;\r
+                   }\r
+\r
+                   ilast = indx[0];\r
+\r
+\r
+                 } else {\r
+                   /* ilast is still OK */\r
+\r
+                   indx[0] = ilast;\r
+\r
+                 }\r
+\r
+                 return(CV_SUCCESS);\r
+\r
+\r
+               }\r
+\r
+               /*\r
+                * CVApolynomialGetY ( -> IMget )\r
+                *\r
+                * This routine uses polynomial interpolation for the forward solution vector. \r
+                * It is typically called by the wrapper routines before calling\r
+                * user provided routines (fB, djacB, bjacB, jtimesB, psolB)) but\r
+                * can be directly called by the user through CVodeGetAdjY.\r
+                */\r
+\r
+               protected int CVApolynomialGetY(CVodeMemRec cv_mem, double t,\r
+                                            NVector y, NVector yS[])\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 DtpntMemRec dt_mem[];\r
+                 PolynomialDataMemRec content;\r
+\r
+                 int flag, dir, order, i, j, is, NS;\r
+                 long indx[] = new long[1]; \r
+                 long base;\r
+                 boolean newpoint[] = new boolean[1];\r
+                 double dt, factor;\r
+\r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+                 dt_mem = ca_mem.dt_mem;\r
+                 \r
+                 /* Local value of Ns */\r
+                \r
+                 NS = (ca_mem.ca_IMinterpSensi && (yS != null)) ? cv_mem.cv_Ns : 0;\r
+\r
+                 /* Get the index in dt_mem */\r
+\r
+                 flag = CVAfindIndex(cv_mem, t, indx, newpoint);\r
+                 if (flag != CV_SUCCESS) return(flag);\r
+\r
+                 /* If we are beyond the left limit but close enough,\r
+                    then return y at the left limit. */\r
+\r
+                 if (indx[0] == 0) {\r
+                   content = dt_mem[0].polynomialContent;\r
+                   N_VScale_Serial(ONE, content.y, y);\r
+                   for (is=0; is<NS; is++) N_VScale_Serial(ONE, content.yS[is], yS[is]);\r
+                   return(CV_SUCCESS);\r
+                 }\r
+\r
+                 /* Scaling factor */\r
+\r
+                 dt = Math.abs(dt_mem[(int)indx[0]].t - dt_mem[(int)(indx[0]-1)].t);\r
+\r
+                 /* Find the direction of the forward integration */\r
+\r
+                 dir = (ca_mem.ca_tfinal - ca_mem.ca_tinitial > ZERO) ? 1 : -1;\r
+\r
+                 /* Establish the base point depending on the integration direction.\r
+                    Modify the base if there are not enough points for the current order */\r
+\r
+                 if (dir == 1) {\r
+                   base = indx[0];\r
+                   content = dt_mem[(int)base].polynomialContent;\r
+                   order = content.order;\r
+                   if(indx[0] < order) base += order-indx[0];\r
+                 } else {\r
+                   base = indx[0]-1;\r
+                   content = dt_mem[(int)base].polynomialContent;\r
+                   order = content.order;\r
+                   if (ca_mem.ca_np-indx[0] > order) base -= indx[0]+order-ca_mem.ca_np;\r
+                 }\r
+\r
+                 /* Recompute Y (divided differences for Newton polynomial) if needed */\r
+\r
+                 if (newpoint[0]) {\r
+\r
+                   /* Store 0-th order DD */\r
+                   if (dir == 1) {\r
+                     for(j=0;j<=order;j++) {\r
+                       ca_mem.ca_T[j] = dt_mem[(int)(base-j)].t;\r
+                       content = dt_mem[(int)(base-j)].polynomialContent;\r
+                       N_VScale_Serial(ONE, content.y, ca_mem.ca_Y[j]);\r
+                       for (is=0; is<NS; is++) N_VScale_Serial(ONE, content.yS[is], ca_mem.ca_YS[j][is]);\r
+                     }\r
+                   } else {\r
+                     for(j=0;j<=order;j++) {\r
+                       ca_mem.ca_T[j] = dt_mem[(int)(base-1+j)].t;\r
+                       content = dt_mem[(int)(base-1+j)].polynomialContent;\r
+                       N_VScale_Serial(ONE, content.y, ca_mem.ca_Y[j]);\r
+                       for (is=0; is<NS; is++) N_VScale_Serial(ONE, content.yS[is], ca_mem.ca_YS[j][is]);\r
+                     }\r
+                   }\r
+\r
+                   /* Compute higher-order DD */\r
+                   for(i=1;i<=order;i++) {\r
+                     for(j=order;j>=i;j--) {\r
+                       factor = dt/(ca_mem.ca_T[j]-ca_mem.ca_T[j-i]);\r
+                       N_VLinearSum_Serial(factor, ca_mem.ca_Y[j], -factor, ca_mem.ca_Y[j-1], ca_mem.ca_Y[j]);\r
+                       for (is=0; is<NS; is++) \r
+                               N_VLinearSum_Serial(factor, ca_mem.ca_YS[j][is], -factor, ca_mem.ca_YS[j-1][is], ca_mem.ca_YS[j][is]);\r
+                     }\r
+                   }\r
+                 }\r
+\r
+                 /* Perform the actual interpolation using nested multiplications */\r
+\r
+                 N_VScale_Serial(ONE, ca_mem.ca_Y[order], y);\r
+                 for (is=0; is<NS; is++) N_VScale_Serial(ONE, ca_mem.ca_YS[order][is], yS[is]);\r
+                 for (i=order-1; i>=0; i--) {\r
+                   factor = (t-ca_mem.ca_T[i])/dt;\r
+                   N_VLinearSum_Serial(factor, y, ONE, ca_mem.ca_Y[i], y);\r
+                   for (is=0; is<NS; is++) N_VLinearSum_Serial(factor, yS[is], ONE, ca_mem.ca_YS[i][is], yS[is]);\r
+                 }\r
+\r
+                 return(CV_SUCCESS);\r
+\r
+               }\r
+               \r
+               /*\r
+                * CVArhs\r
+                *\r
+                * This routine interfaces to the CVRhsFnB (or CVRhsFnBS) routine \r
+                * provided by the user.\r
+                */\r
+\r
+\r
+               private int CVArhs(double t, NVector yB, \r
+                                 NVector yBdot, CVodeMemRec cv_mem)\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 CVodeBMemRec cvB_mem;\r
+                 int flag = 0;\r
+                 int retval = 0;\r
+\r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                 cvB_mem = ca_mem.ca_bckpbCrt;\r
+\r
+                 /* Get forward solution from interpolation */\r
+\r
+                 if (ca_mem.ca_IMinterpSensi)\r
+                       if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                       flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, ca_mem.ca_yStmp);\r
+                       }\r
+                       else {\r
+                               flag = CVApolynomialGetY(cv_mem, t, ca_mem.ca_ytmp, ca_mem.ca_yStmp);   \r
+                       }\r
+                 else \r
+                         if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                         flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                         }\r
+                         else {\r
+                                 flag = CVApolynomialGetY(cv_mem, t, ca_mem.ca_ytmp, null);  \r
+                         }\r
+\r
+                 if (flag != CV_SUCCESS) {\r
+                   cvProcessError(cv_mem, -1, "CVODEA", "CVArhs", MSGCV_BAD_TINTERP, t);\r
+                   return(-1);\r
+                 }\r
+\r
+                 /* Call the user's RHS function */\r
+\r
+                 if (cvB_mem.cv_f_withSensi)\r
+                       if (testMode) {\r
+                               //retval = fBS1TestMode(t, ca_mem.ca_ytmp, ca_mem.ca_yStmp, yB, yBdot, cvB_mem.cv_user_data);   \r
+                       }\r
+                       else {\r
+                       //retval = fBS1(t, ca_mem.ca_ytmp, ca_mem.ca_yStmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                       }\r
+                 else\r
+                         if (testMode) {\r
+                                 retval = fBTestMode(t, ca_mem.ca_ytmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                         }\r
+                         else {\r
+                         retval = fB(t, ca_mem.ca_ytmp, yB, yBdot, cvB_mem.cv_user_data);\r
+                         }\r
+\r
+                 return(retval);\r
+               }\r
+               \r
+               /*\r
+                * CVArhsQ\r
+                *\r
+                * This routine interfaces to the CVQuadRhsFnB (or CVQuadRhsFnBS) routine\r
+                * provided by the user.\r
+                */\r
+\r
+       private int CVArhsQ(double t, NVector yB, \r
+                                  NVector qBdot, CVodeMemRec cv_mem)\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 CVodeBMemRec cvB_mem;\r
+                 /* int flag; */\r
+                 int retval = 0;\r
+\r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                 cvB_mem = ca_mem.ca_bckpbCrt;\r
+\r
+                 /* Get forward solution from interpolation */\r
+\r
+                 if (ca_mem.ca_IMinterpSensi) {\r
+                         if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
+                          /* flag = */ CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, ca_mem.ca_yStmp);\r
+                         }\r
+                         else {\r
+                                 /* flag = */ CVApolynomialGetY(cv_mem, t, ca_mem.ca_ytmp, ca_mem.ca_yStmp);\r
+                         }\r
+                 }\r
+                 else {\r
+                         if (ca_mem.ca_IMget == CVAhermiteGetY_select)  {\r
+                           /* flag = */ CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                         }\r
+                         else {\r
+                                 /* flag = */ CVApolynomialGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
+                         }\r
+                 }\r
+\r
+                 /* Call the user's RHS function */\r
+\r
+                 if (cvB_mem.cv_fQ_withSensi) {\r
+                   //retval = (cvB_mem->cv_fQs)(t, ca_mem->ca_ytmp, ca_mem->ca_yStmp, yB, qBdot, cvB_mem->cv_user_data);\r
+                 }\r
+                 else\r
+                       if (testMode) {\r
+                               retval = fQBTestMode(t, ca_mem.ca_ytmp, yB, qBdot, cvB_mem.cv_user_data);\r
+                       }\r
+                       else {\r
+                   retval = fQB(t, ca_mem.ca_ytmp, yB, qBdot, cvB_mem.cv_user_data);\r
+                       }\r
+\r
+                 return(retval);\r
+               }\r
+\r
 \r
 }
\ No newline at end of file