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

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

index babd29a..a2f8050 100644 (file)
@@ -4,7 +4,7 @@ package gov.nih.mipav.model.algorithms;
 import gov.nih.mipav.view.MipavUtil;\r
 \r
 public abstract class CVODES_ASA extends CVODES {\r
-       final int CVArhs_select = 100;\r
+       \r
        public CVODES_ASA() {\r
                super();\r
                runcvsRoberts_ASAi_dns();\r
@@ -141,8 +141,10 @@ public abstract class CVODES_ASA extends CVODES {
                int flag;\r
         int f = cvsRoberts_ASAi_dns;\r
         int Jac = cvsRoberts_ASAi_dns;\r
+        int JacB = cvsRoberts_ASAi_dns;\r
                int ewt_select = cvEwtUser_select1;\r
                int fQ= cvsRoberts_ASAi_dns;\r
+               int fQB = cvsRoberts_ASAi_dns;\r
                double A[][];\r
                SUNLinearSolver LS;\r
                int steps;\r
@@ -156,6 +158,8 @@ public abstract class CVODES_ASA extends CVODES {
                double abstolQB;\r
                int indexB[] = new int[1];\r
                int fB = cvsRoberts_ASAi_dns;\r
+               double AB[][];\r
+               SUNLinearSolver LSB;\r
                \r
                // Print problem description \r
                System.out.printf("\nAdjoint Sensitivity Example for Chemical Kinetics\n");\r
@@ -358,7 +362,74 @@ public abstract class CVODES_ASA extends CVODES {
                // Call CVodeInitB to allocate internal memory and initialize the\r
                // backward problem.\r
                flag = CVodeInitB(cvode_mem, indexB[0], fB, TB1, yB);\r
-       }\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the scalar relative and absolute tolerances.\r
+               flag = CVodeSStolerancesB(cvode_mem, indexB[0], reltolB, abstolB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               // Attach the user data for backward problem\r
+               flag = CVodeSetUserDataB(cvode_mem, indexB[0], data);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Create dense SUNMatrix for user in linear solver\r
+               AB = new double[NEQ][NEQ];\r
+               \r
+               // Create dense SUNLinearSolver object\r
+               LSB = SUNDenseLinearSolver(yB, AB);\r
+               if (LSB == null) {\r
+                       return;\r
+               }\r
+               \r
+               // Attach the matrix and linear solver\r
+               flag = CVDlsSetLinearSolverB(cvode_mem, indexB[0], LSB, AB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the user-supplied Jacobian routine JacB\r
+               flag = CVDlsSetJacFnB(cvode_mem, indexB[0], JacB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeQuadInitB to allocate internal memory and initialize backward\r
+               // quadrature integration.\r
+               flag = CVodeQuadInitB(cvode_mem, indexB[0], fQB, qB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeSetQuadErrConB to specify whether or not the quadrature variables\r
+               // are to be used in the step size control mechanism within CVODES.  Call\r
+               // CVodeQuadSStoerancesB or CVodeQuadSVtolerancesB to specify the integration\r
+               // tolerances for the quadrature variables.\r
+               flag = CVodeSetQuadErrConB(cvode_mem, indexB[0], true);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeQUadSStolerancesB to specify the scalar relative and absolute tolerances \r
+               // for the backward problem.\r
+               flag = CVodeQuadSStolerancesB(cvode_mem, indexB[0], reltolB, abstolQB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Backward integration\r
+               System.out.printf("Backward integration from tB0 = " + TB1 + "\n\n");\r
+               \r
+               // First get results at t = TBout1\r
+               \r
+               // Call CVodeBto integrate the backward ODE problem.\r
+               flag = CVodeB(cvode_mem, TBout1, CV_NORMAL);\r
+       } \r
        \r
        /*\r
         * CVodeQuadInit\r
@@ -679,1392 +750,1836 @@ public abstract class CVODES_ASA extends CVODES {
        } \r
 \r
        \r
-               /*\r
-                * CVodeF\r
-                *\r
-                * This routine integrates to tout and returns solution into yout.\r
-                * In the same time, it stores check point data every 'steps' steps. \r
-                * \r
-                * CVodeF can be called repeatedly by the user.\r
-                *\r
-                * ncheckPtr points to the number of check points stored so far.\r
-                */\r
-\r
-               private int CVodeF(CVodeMemRec cv_mem, double tout, NVector yout, \r
-                          double tret[], int itask, int ncheckPtr[])\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 CkpntMemRec tmp;\r
-                 DtpntMemRec dt_mem[];\r
-                 int flag, i;\r
-                 boolean iret, allocOK;\r
-\r
-                 /* Check if cvode_mem exists */\r
-                 if (cv_mem == null) {\r
-                   cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeF", MSGCV_NO_MEM);\r
-                   return(CV_MEM_NULL);\r
-                 }\r
-\r
-                 /* Was ASA initialized? */\r
-                 if (cv_mem.cv_adjMallocDone == false) {\r
-                   cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeF", MSGCV_NO_ADJ);\r
-                   return(CV_NO_ADJ);\r
-                 } \r
-\r
-                 ca_mem = cv_mem.cv_adj_mem;\r
-\r
-                 /* Check for yout != NULL */\r
-                 if (yout == null) {\r
-                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_YOUT_NULL);\r
-                   return(CV_ILL_INPUT);\r
-                 }\r
-                 \r
-                 /* Check for tret != NULL */\r
-                 if (tret == null) {\r
-                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_TRET_NULL);\r
-                   return(CV_ILL_INPUT);\r
-                 }\r
-\r
-                 /* Check for valid itask */\r
-                 if ( (itask != CV_NORMAL) && (itask != CV_ONE_STEP) ) {\r
-                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_BAD_ITASK);\r
-                   return(CV_ILL_INPUT);\r
-                 }\r
-\r
-                 /* All error checking done */\r
-\r
-                 dt_mem = ca_mem.dt_mem;\r
+       /*\r
+        * CVodeF\r
+        *\r
+        * This routine integrates to tout and returns solution into yout.\r
+        * In the same time, it stores check point data every 'steps' steps. \r
+        * \r
+        * CVodeF can be called repeatedly by the user.\r
+        *\r
+        * ncheckPtr points to the number of check points stored so far.\r
+        */\r
 \r
-                 /* If tstop is enabled, store some info */\r
-                 if (cv_mem.cv_tstopset) {\r
-                   ca_mem.ca_tstopCVodeFcall = true;\r
-                   ca_mem.ca_tstopCVodeF = cv_mem.cv_tstop;\r
-                 }\r
+       private int CVodeF(CVodeMemRec cv_mem, double tout, NVector yout, \r
+                  double tret[], int itask, int ncheckPtr[])\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CkpntMemRec tmp;\r
+         DtpntMemRec dt_mem[];\r
+         int flag, i;\r
+         boolean iret, allocOK;\r
 \r
-                 /* We will call CVode in CV_ONE_STEP mode, regardless\r
-                  * of what itask is, so flag if we need to return */\r
-                 if (itask == CV_ONE_STEP) iret = true;\r
-                 else                      iret = false;\r
-\r
-                 /* On the first step:\r
-                  *   - set tinitial\r
-                  *   - initialize list of check points\r
-                  *   - if needed, initialize the interpolation module\r
-                  *   - load dt_mem[0]\r
-                  * On subsequent steps, test if taking a new step is necessary. \r
-                  */\r
-                 if ( ca_mem.ca_firstCVodeFcall ) {\r
-\r
-                   ca_mem.ca_tinitial = cv_mem.cv_tn;\r
-\r
-                   ca_mem.ck_mem = CVAckpntInit(cv_mem);\r
-                   if (ca_mem.ck_mem == null) {\r
-                     cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
-                     return(CV_MEM_FAIL);\r
-                   }\r
-\r
-                   if ( !ca_mem.ca_IMmallocDone ) {\r
-\r
-                     /* Do we need to store sensitivities? */\r
-                     if (!cv_mem.cv_sensi) ca_mem.ca_IMstoreSensi = false;\r
-\r
-                     /* Allocate space for interpolation data */\r
-                     if (ca_mem.ca_IMmalloc == CVAhermiteMalloc_select) {\r
-                         allocOK = CVAhermiteMalloc(cv_mem);\r
-                     }\r
-                     else {\r
-                         allocOK = CVApolynomialMalloc(cv_mem);\r
-                     }\r
-                     if (!allocOK) {\r
-                       cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
-                       return(CV_MEM_FAIL);\r
-                     }\r
-\r
-                     /* Rename zn and, if needed, znS for use in interpolation */\r
-                     for (i=0;i<L_MAX;i++) ca_mem.ca_Y[i] = cv_mem.cv_zn[i];\r
-                     if (ca_mem.ca_IMstoreSensi) {\r
-                       for (i=0;i<L_MAX;i++) ca_mem.ca_YS[i] = cv_mem.cv_znS[i];\r
-                     }\r
-\r
-                     ca_mem.ca_IMmallocDone = true;\r
-\r
-                   }\r
-\r
-                   dt_mem[0].t = ca_mem.ck_mem.ck_t0;\r
-                   if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
-                       CVAhermiteStorePnt(cv_mem, dt_mem[0]);\r
-                   }\r
-                   else {\r
-                       CVApolynomialStorePnt(cv_mem, dt_mem[0]);       \r
-                   }\r
-\r
-                   ca_mem.ca_firstCVodeFcall = false;\r
-\r
-                 } else if ( (cv_mem.cv_tn - tout)*cv_mem.cv_h >= ZERO ) {\r
-\r
-                   /* If tout was passed, return interpolated solution. \r
-                      No changes to ck_mem or dt_mem are needed. */\r
-                   tret[0] = tout;\r
-                   flag = CVodeGetDky(cv_mem, tout, 0, yout);\r
-                   ncheckPtr[0] = ca_mem.ca_nckpnts;\r
-                   ca_mem.ca_IMnewData = true;\r
-                   ca_mem.ca_ckpntData = ca_mem.ck_mem;\r
-                   ca_mem.ca_np = cv_mem.cv_nst % ca_mem.ca_nsteps + 1;\r
-\r
-                   return(flag);\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeF", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
 \r
-                 }\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeF", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
 \r
-                 /* Integrate to tout (in CV_ONE_STEP mode) while loading check points */\r
-                 for(;;) {\r
-\r
-                   /* Perform one step of the integration */\r
-\r
-                   flag = CVode(cv_mem, tout, yout, tret, CV_ONE_STEP);\r
-                   if (flag < 0) break;\r
-\r
-                   /* Test if a new check point is needed */\r
-\r
-                   if ( cv_mem.cv_nst % ca_mem.ca_nsteps == 0 ) {\r
-\r
-                     ca_mem.ck_mem.ck_t1 = tret[0];\r
-\r
-                     /* Create a new check point, load it, and append it to the list */\r
-                     tmp = CVAckpntNew(cv_mem);\r
-                     if (tmp == null) {\r
-                       cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
-                       flag = CV_MEM_FAIL;\r
-                       break;\r
-                     }\r
-                     tmp.ck_next = ca_mem.ck_mem;\r
-                     ca_mem.ck_mem = tmp;\r
-                     ca_mem.ca_nckpnts++;\r
-                     cv_mem.cv_forceSetup = true;\r
-                     \r
-                     /* Reset i=0 and load dt_mem[0] */\r
-                     dt_mem[0].t = ca_mem.ck_mem.ck_t0;\r
-                     if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
-                         CVAhermiteStorePnt(cv_mem, dt_mem[0]);\r
-                     }\r
-                     else {\r
-                         CVApolynomialStorePnt(cv_mem, dt_mem[0]);     \r
-                     }\r
-\r
-                   } else {\r
-\r
-                     /* Load next point in dt_mem */\r
-                     dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)].t = tret[0];\r
-                     if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
-                         CVAhermiteStorePnt(cv_mem, dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)]);\r
-                     }\r
-                     else {\r
-                         CVApolynomialStorePnt(cv_mem, dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)]);  \r
-                     }\r
-\r
-                   }\r
-\r
-                   /* Set t1 field of the current ckeck point structure\r
-                      for the case in which there will be no future\r
-                      check points */\r
-                   ca_mem.ck_mem.ck_t1 = tret[0];\r
-\r
-                   /* tfinal is now set to *tret */\r
-                   ca_mem.ca_tfinal = tret[0];\r
-\r
-                   /* Return if in CV_ONE_STEP mode */\r
-                   if (iret) break;\r
-\r
-                   /* Return if tout reached */\r
-                   if ( (tret[0] - tout)*cv_mem.cv_h >= ZERO ) {\r
-                     tret[0] = tout;\r
-                     CVodeGetDky(cv_mem, tout, 0, yout);\r
-                     /* Reset tretlast in cv_mem so that CVodeGetQuad and CVodeGetSens \r
-                      * evaluate quadratures and/or sensitivities at the proper time */\r
-                     cv_mem.cv_tretlast = tout;\r
-                     break;\r
-                   }\r
-\r
-                 } /* end of for(;;)() */\r
-\r
-                 /* Get ncheck from ca_mem */ \r
-                 ncheckPtr[0] = ca_mem.ca_nckpnts;\r
-\r
-                 /* Data is available for the last interval */\r
-                 ca_mem.ca_IMnewData = true;\r
-                 ca_mem.ca_ckpntData = ca_mem.ck_mem;\r
-                 ca_mem.ca_np = cv_mem.cv_nst % ca_mem.ca_nsteps + 1;\r
-\r
-                 return(flag);\r
-               }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-               /*\r
-                * CVAckpntInit\r
-                *\r
-                * This routine initializes the check point linked list with \r
-                * information from the initial time.\r
-                */\r
-\r
-               private CkpntMemRec CVAckpntInit(CVodeMemRec cv_mem)\r
-               {\r
-                 CkpntMemRec ck_mem;\r
-                 int is;\r
-\r
-                 /* Allocate space for ckdata */\r
-                 ck_mem = null;\r
-                 ck_mem = new CkpntMemRec();\r
-\r
-                 ck_mem.ck_zn[0] = N_VClone_Serial(cv_mem.cv_tempv);\r
-                 if (ck_mem.ck_zn[0] == null) {\r
-                   ck_mem = null;\r
-                   return(null);\r
-                 }\r
-                 \r
-                 ck_mem.ck_zn[1] = N_VClone_Serial(cv_mem.cv_tempv);\r
-                 if (ck_mem.ck_zn[1] == null) {\r
-                   N_VDestroy(ck_mem.ck_zn[0]);\r
-                   ck_mem = null;\r
-                   return(null);\r
-                 }\r
+         /* Check for yout != NULL */\r
+         if (yout == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_YOUT_NULL);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         \r
+         /* Check for tret != NULL */\r
+         if (tret == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_TRET_NULL);\r
+           return(CV_ILL_INPUT);\r
+         }\r
 \r
-                 /* ck_mem->ck_zn[qmax] was not allocated */\r
-                 ck_mem.ck_zqm = 0;\r
+         /* Check for valid itask */\r
+         if ( (itask != CV_NORMAL) && (itask != CV_ONE_STEP) ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_BAD_ITASK);\r
+           return(CV_ILL_INPUT);\r
+         }\r
 \r
-                 /* Load ckdata from cv_mem */\r
-                 N_VScale_Serial(ONE, cv_mem.cv_zn[0], ck_mem.ck_zn[0]);\r
-                 ck_mem.ck_t0    = cv_mem.cv_tn;\r
-                 ck_mem.ck_nst   = 0;\r
-                 ck_mem.ck_q     = 1;\r
-                 ck_mem.ck_h     = 0.0;\r
-                 \r
-                 /* Do we need to carry quadratures */\r
-                 ck_mem.ck_quadr = cv_mem.cv_quadr && cv_mem.cv_errconQ;\r
+         /* All error checking done */\r
 \r
-                 if (ck_mem.ck_quadr) {\r
+         dt_mem = ca_mem.dt_mem;\r
 \r
-                   ck_mem.ck_znQ[0] = N_VClone_Serial(cv_mem.cv_tempvQ);\r
-                   if (ck_mem.ck_znQ[0] == null) {\r
-                     N_VDestroy(ck_mem.ck_zn[0]);\r
-                     N_VDestroy(ck_mem.ck_zn[1]);\r
-                     ck_mem = null;\r
-                     return(null);\r
-                   }\r
+         /* If tstop is enabled, store some info */\r
+         if (cv_mem.cv_tstopset) {\r
+           ca_mem.ca_tstopCVodeFcall = true;\r
+           ca_mem.ca_tstopCVodeF = cv_mem.cv_tstop;\r
+         }\r
 \r
-                   N_VScale_Serial(ONE, cv_mem.cv_znQ[0], ck_mem.ck_znQ[0]);\r
+         /* We will call CVode in CV_ONE_STEP mode, regardless\r
+          * of what itask is, so flag if we need to return */\r
+         if (itask == CV_ONE_STEP) iret = true;\r
+         else                      iret = false;\r
+\r
+         /* On the first step:\r
+          *   - set tinitial\r
+          *   - initialize list of check points\r
+          *   - if needed, initialize the interpolation module\r
+          *   - load dt_mem[0]\r
+          * On subsequent steps, test if taking a new step is necessary. \r
+          */\r
+         if ( ca_mem.ca_firstCVodeFcall ) {\r
+\r
+           ca_mem.ca_tinitial = cv_mem.cv_tn;\r
+\r
+           ca_mem.ck_mem = CVAckpntInit(cv_mem);\r
+           if (ca_mem.ck_mem == null) {\r
+             cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
+             return(CV_MEM_FAIL);\r
+           }\r
 \r
-                 }\r
+           if ( !ca_mem.ca_IMmallocDone ) {\r
 \r
-                 /* Do we need to carry sensitivities? */\r
-                 ck_mem.ck_sensi = cv_mem.cv_sensi;\r
+             /* Do we need to store sensitivities? */\r
+             if (!cv_mem.cv_sensi) ca_mem.ca_IMstoreSensi = false;\r
 \r
-                 if (ck_mem.ck_sensi) {\r
+             /* Allocate space for interpolation data */\r
+             if (ca_mem.ca_IMmalloc == CVAhermiteMalloc_select) {\r
+                 allocOK = CVAhermiteMalloc(cv_mem);\r
+             }\r
+             else {\r
+                 allocOK = CVApolynomialMalloc(cv_mem);\r
+             }\r
+             if (!allocOK) {\r
+               cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
+               return(CV_MEM_FAIL);\r
+             }\r
 \r
-                   ck_mem.ck_Ns = cv_mem.cv_Ns;\r
+             /* Rename zn and, if needed, znS for use in interpolation */\r
+             for (i=0;i<L_MAX;i++) ca_mem.ca_Y[i] = cv_mem.cv_zn[i];\r
+             if (ca_mem.ca_IMstoreSensi) {\r
+               for (i=0;i<L_MAX;i++) ca_mem.ca_YS[i] = cv_mem.cv_znS[i];\r
+             }\r
 \r
-                   ck_mem.ck_znS[0] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                   if (ck_mem.ck_znS[0] == null) {\r
-                     N_VDestroy(ck_mem.ck_zn[0]);\r
-                     N_VDestroy(ck_mem.ck_zn[1]);\r
-                     if (ck_mem.ck_quadr) N_VDestroy(ck_mem.ck_znQ[0]);\r
-                     ck_mem = null;\r
-                     return(null);\r
-                   }\r
+             ca_mem.ca_IMmallocDone = true;\r
 \r
-                   for (is=0; is<cv_mem.cv_Ns; is++)\r
-                     N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], ck_mem.ck_znS[0][is]);\r
+           }\r
 \r
-                 }\r
+           dt_mem[0].t = ca_mem.ck_mem.ck_t0;\r
+           if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
+               CVAhermiteStorePnt(cv_mem, dt_mem[0]);\r
+           }\r
+           else {\r
+               CVApolynomialStorePnt(cv_mem, dt_mem[0]);       \r
+           }\r
 \r
-                 /* Do we need to carry quadrature sensitivities? */\r
-                 ck_mem.ck_quadr_sensi = cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS;\r
-\r
-                 if (ck_mem.ck_quadr_sensi) {\r
-                   ck_mem.ck_znQS[0] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
-                   if (ck_mem.ck_znQS[0] == null) {\r
-                     N_VDestroy(ck_mem.ck_zn[0]);\r
-                     N_VDestroy(ck_mem.ck_zn[1]);\r
-                     if (ck_mem.ck_quadr) N_VDestroy(ck_mem.ck_znQ[0]);\r
-                     N_VDestroyVectorArray_Serial(ck_mem.ck_znS[0], cv_mem.cv_Ns);\r
-                     ck_mem = null;\r
-                     return(null);\r
-                   }\r
-                   \r
-                   for (is=0; is<cv_mem.cv_Ns; is++)\r
-                     N_VScale_Serial(ONE, cv_mem.cv_znQS[0][is], ck_mem.ck_znQS[0][is]);\r
+           ca_mem.ca_firstCVodeFcall = false;\r
 \r
-                 }\r
+         } else if ( (cv_mem.cv_tn - tout)*cv_mem.cv_h >= ZERO ) {\r
 \r
-                 /* Next in list */\r
-                 ck_mem.ck_next  = null;\r
+           /* If tout was passed, return interpolated solution. \r
+              No changes to ck_mem or dt_mem are needed. */\r
+           tret[0] = tout;\r
+           flag = CVodeGetDky(cv_mem, tout, 0, yout);\r
+           ncheckPtr[0] = ca_mem.ca_nckpnts;\r
+           ca_mem.ca_IMnewData = true;\r
+           ca_mem.ca_ckpntData = ca_mem.ck_mem;\r
+           ca_mem.ca_np = cv_mem.cv_nst % ca_mem.ca_nsteps + 1;\r
 \r
-                 return(ck_mem);\r
-               }\r
-               \r
-               /*\r
-                * CVAhermiteMalloc\r
-                *\r
-                * This routine allocates memory for storing information at all\r
-                * intermediate points between two consecutive check points. \r
-                * This data is then used to interpolate the forward solution \r
-                * at any other time.\r
-                */\r
-\r
-               private boolean CVAhermiteMalloc(CVodeMemRec cv_mem)\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 DtpntMemRec dt_mem[];\r
-                 HermiteDataMemRec content;\r
-                 int i, ii=0;\r
-                 boolean allocOK;\r
-\r
-                 allocOK = true;\r
-\r
-                 ca_mem = cv_mem.cv_adj_mem;\r
-\r
-                 /* Allocate space for the vectors ytmp and yStmp */\r
-\r
-                 ca_mem.ca_ytmp = N_VClone(cv_mem.cv_tempv);\r
-                 if (ca_mem.ca_ytmp == null) {\r
-                   return(false);\r
-                 }\r
+           return(flag);\r
 \r
-                 if (ca_mem.ca_IMstoreSensi) {\r
-                   ca_mem.ca_yStmp = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                   if (ca_mem.ca_yStmp == null) {\r
-                     N_VDestroy(ca_mem.ca_ytmp);\r
-                     return(false);\r
-                   }\r
-                 }\r
+         }\r
 \r
-                 /* Allocate space for the content field of the dt structures */\r
-\r
-                 dt_mem = ca_mem.dt_mem;\r
-\r
-                 for (i=0; i<=ca_mem.ca_nsteps; i++) {\r
-\r
-                   content = null;\r
-                   content = new HermiteDataMemRec();\r
-\r
-                   content.y = N_VClone(cv_mem.cv_tempv);\r
-                   if (content.y == null) {\r
-                     content = null;\r
-                     ii = i;\r
-                     allocOK = false;\r
-                     break;\r
-                   }\r
-\r
-                   content.yd = N_VClone(cv_mem.cv_tempv);\r
-                   if (content.yd == null) {\r
-                     N_VDestroy(content.y);\r
-                     content = null;\r
-                     ii = i;\r
-                     allocOK = false;\r
-                     break;\r
-                   }\r
-\r
-                   if (ca_mem.ca_IMstoreSensi) {\r
-\r
-                     content.yS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                     if (content.yS == null) {\r
-                       N_VDestroy(content.y);\r
-                       N_VDestroy(content.yd);\r
-                       content = null;\r
-                       ii = i;\r
-                       allocOK = false;\r
-                       break;\r
-                     }\r
-\r
-                     content.ySd = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                     if (content.ySd == null) {\r
-                       N_VDestroy(content.y);\r
-                       N_VDestroy(content.yd);\r
-                       N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
-                       content = null;\r
-                       ii = i;\r
-                       allocOK = false;\r
-                       break;\r
-                     }\r
-                     \r
-                   }\r
-                   \r
-                   dt_mem[i].content = CV_HERMITE;\r
-                   dt_mem[i].hermiteContent = content;\r
-\r
-                 } \r
-\r
-                 /* If an error occurred, deallocate and return */\r
-\r
-                 if (!allocOK) {\r
-\r
-                   N_VDestroy(ca_mem.ca_ytmp);\r
-\r
-                   if (ca_mem.ca_IMstoreSensi) {\r
-                     N_VDestroyVectorArray_Serial(ca_mem.ca_yStmp, cv_mem.cv_Ns);\r
-                   }\r
-\r
-                   for (i=0; i<ii; i++) {\r
-                     content = (dt_mem[i].hermiteContent);\r
-                     N_VDestroy(content.y);\r
-                     N_VDestroy(content.yd);\r
-                     if (ca_mem.ca_IMstoreSensi) {\r
-                       N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
-                       N_VDestroyVectorArray_Serial(content.ySd, cv_mem.cv_Ns);\r
-                     }\r
-                     dt_mem[i].hermiteContent = null;\r
-                   }\r
+         /* Integrate to tout (in CV_ONE_STEP mode) while loading check points */\r
+         for(;;) {\r
 \r
-                 }\r
+           /* Perform one step of the integration */\r
 \r
-                 return(allocOK);\r
-               }\r
+           flag = CVode(cv_mem, tout, yout, tret, CV_ONE_STEP);\r
+           if (flag < 0) break;\r
 \r
-               \r
-               /*\r
-                * CVApolynomialMalloc\r
-                *\r
-                * This routine allocates memory for storing information at all\r
-                * intermediate points between two consecutive check points. \r
-                * This data is then used to interpolate the forward solution \r
-                * at any other time.\r
-                */\r
-\r
-               private boolean CVApolynomialMalloc(CVodeMemRec cv_mem)\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 DtpntMemRec dt_mem[];\r
-                 PolynomialDataMemRec content;\r
-                 int i, ii=0;\r
-                 boolean allocOK;\r
-\r
-                 allocOK = true;\r
-\r
-                 ca_mem = cv_mem.cv_adj_mem;\r
-\r
-                 /* Allocate space for the vectors ytmp and yStmp */\r
-\r
-                 ca_mem.ca_ytmp = N_VClone(cv_mem.cv_tempv);\r
-                 if (ca_mem.ca_ytmp == null) {\r
-                   return(false);\r
-                 }\r
+           /* Test if a new check point is needed */\r
 \r
-                 if (ca_mem.ca_IMstoreSensi) {\r
-                   ca_mem.ca_yStmp = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                   if (ca_mem.ca_yStmp == null) {\r
-                     N_VDestroy(ca_mem.ca_ytmp);\r
-                     return(false);\r
-                   }\r
-                 }\r
+           if ( cv_mem.cv_nst % ca_mem.ca_nsteps == 0 ) {\r
 \r
-                 /* Allocate space for the content field of the dt structures */\r
+             ca_mem.ck_mem.ck_t1 = tret[0];\r
 \r
-                 dt_mem = ca_mem.dt_mem;\r
+             /* Create a new check point, load it, and append it to the list */\r
+             tmp = CVAckpntNew(cv_mem);\r
+             if (tmp == null) {\r
+               cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);\r
+               flag = CV_MEM_FAIL;\r
+               break;\r
+             }\r
+             tmp.ck_next = ca_mem.ck_mem;\r
+             ca_mem.ck_mem = tmp;\r
+             ca_mem.ca_nckpnts++;\r
+             cv_mem.cv_forceSetup = true;\r
+             \r
+             /* Reset i=0 and load dt_mem[0] */\r
+             dt_mem[0].t = ca_mem.ck_mem.ck_t0;\r
+             if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
+                 CVAhermiteStorePnt(cv_mem, dt_mem[0]);\r
+             }\r
+             else {\r
+                 CVApolynomialStorePnt(cv_mem, dt_mem[0]);     \r
+             }\r
 \r
-                 for (i=0; i<=ca_mem.ca_nsteps; i++) {\r
+           } else {\r
 \r
-                   content = null;\r
-                   content = new PolynomialDataMemRec();\r
+             /* Load next point in dt_mem */\r
+             dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)].t = tret[0];\r
+             if (ca_mem.ca_IMstore == CVAhermiteStorePnt_select) {\r
+                 CVAhermiteStorePnt(cv_mem, dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)]);\r
+             }\r
+             else {\r
+                 CVApolynomialStorePnt(cv_mem, dt_mem[(int)(cv_mem.cv_nst % ca_mem.ca_nsteps)]);  \r
+             }\r
 \r
-                   content.y = N_VClone(cv_mem.cv_tempv);\r
-                   if (content.y == null) {\r
-                     content = null;\r
-                     ii = i;\r
-                     allocOK = false;\r
-                     break;\r
-                   }\r
+           }\r
 \r
-                   if (ca_mem.ca_IMstoreSensi) {\r
+           /* Set t1 field of the current ckeck point structure\r
+              for the case in which there will be no future\r
+              check points */\r
+           ca_mem.ck_mem.ck_t1 = tret[0];\r
+\r
+           /* tfinal is now set to *tret */\r
+           ca_mem.ca_tfinal = tret[0];\r
+\r
+           /* Return if in CV_ONE_STEP mode */\r
+           if (iret) break;\r
+\r
+           /* Return if tout reached */\r
+           if ( (tret[0] - tout)*cv_mem.cv_h >= ZERO ) {\r
+             tret[0] = tout;\r
+             CVodeGetDky(cv_mem, tout, 0, yout);\r
+             /* Reset tretlast in cv_mem so that CVodeGetQuad and CVodeGetSens \r
+              * evaluate quadratures and/or sensitivities at the proper time */\r
+             cv_mem.cv_tretlast = tout;\r
+             break;\r
+           }\r
 \r
-                     content.yS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                     if (content.yS == null) {\r
-                       N_VDestroy(content.y);\r
-                       content = null;\r
-                       ii = i;\r
-                       allocOK = false;\r
-                       break;\r
-                     }\r
+         } /* end of for(;;)() */\r
 \r
-                   }\r
+         /* Get ncheck from ca_mem */ \r
+         ncheckPtr[0] = ca_mem.ca_nckpnts;\r
 \r
-                   dt_mem[i].content = CV_POLYNOMIAL;\r
-                   dt_mem[i].polynomialContent = content;\r
+         /* Data is available for the last interval */\r
+         ca_mem.ca_IMnewData = true;\r
+         ca_mem.ca_ckpntData = ca_mem.ck_mem;\r
+         ca_mem.ca_np = cv_mem.cv_nst % ca_mem.ca_nsteps + 1;\r
 \r
-                 } \r
+         return(flag);\r
+       }\r
 \r
-                 /* If an error occurred, deallocate and return */\r
+       /*\r
+        * CVAckpntInit\r
+        *\r
+        * This routine initializes the check point linked list with \r
+        * information from the initial time.\r
+        */\r
 \r
-                 if (!allocOK) {\r
+       private CkpntMemRec CVAckpntInit(CVodeMemRec cv_mem)\r
+       {\r
+         CkpntMemRec ck_mem;\r
+         int is;\r
 \r
-                   N_VDestroy(ca_mem.ca_ytmp);\r
+         /* Allocate space for ckdata */\r
+         ck_mem = null;\r
+         ck_mem = new CkpntMemRec();\r
 \r
-                   if (ca_mem.ca_IMstoreSensi) {\r
-                     N_VDestroyVectorArray_Serial(ca_mem.ca_yStmp, cv_mem.cv_Ns);\r
-                   }\r
+         ck_mem.ck_zn[0] = N_VClone_Serial(cv_mem.cv_tempv);\r
+         if (ck_mem.ck_zn[0] == null) {\r
+           ck_mem = null;\r
+           return(null);\r
+         }\r
+         \r
+         ck_mem.ck_zn[1] = N_VClone_Serial(cv_mem.cv_tempv);\r
+         if (ck_mem.ck_zn[1] == null) {\r
+           N_VDestroy(ck_mem.ck_zn[0]);\r
+           ck_mem = null;\r
+           return(null);\r
+         }\r
 \r
-                   for (i=0; i<ii; i++) {\r
-                     content = (dt_mem[i].polynomialContent);\r
-                     N_VDestroy(content.y);\r
-                     if (ca_mem.ca_IMstoreSensi) {\r
-                       N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
-                     }\r
-                     dt_mem[i].polynomialContent = null;\r
-                   }\r
+         /* ck_mem->ck_zn[qmax] was not allocated */\r
+         ck_mem.ck_zqm = 0;\r
 \r
-                 }\r
+         /* Load ckdata from cv_mem */\r
+         N_VScale_Serial(ONE, cv_mem.cv_zn[0], ck_mem.ck_zn[0]);\r
+         ck_mem.ck_t0    = cv_mem.cv_tn;\r
+         ck_mem.ck_nst   = 0;\r
+         ck_mem.ck_q     = 1;\r
+         ck_mem.ck_h     = 0.0;\r
+         \r
+         /* Do we need to carry quadratures */\r
+         ck_mem.ck_quadr = cv_mem.cv_quadr && cv_mem.cv_errconQ;\r
 \r
-                 return(allocOK);\r
+         if (ck_mem.ck_quadr) {\r
 \r
-               }\r
-               \r
-               /*\r
-                * CVAhermiteStorePnt ( -> IMstore )\r
-                *\r
-                * This routine stores a new point (y,yd) in the structure d for use\r
-                * in the cubic Hermite interpolation.\r
-                * Note that the time is already stored.\r
-                */\r
-\r
-               private int CVAhermiteStorePnt(CVodeMemRec cv_mem, DtpntMemRec d)\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 HermiteDataMemRec content;\r
-                 int is;\r
-                 /* int retval; */\r
-\r
-                 ca_mem = cv_mem.cv_adj_mem;\r
-\r
-                 content = d.hermiteContent;\r
-\r
-                 /* Load solution */\r
-\r
-                 N_VScale_Serial(ONE, cv_mem.cv_zn[0], content.y);\r
-                 \r
-                 if (ca_mem.ca_IMstoreSensi) {\r
-                   for (is=0; is<cv_mem.cv_Ns; is++) \r
-                     N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], content.yS[is]);\r
-                 }\r
+           ck_mem.ck_znQ[0] = N_VClone_Serial(cv_mem.cv_tempvQ);\r
+           if (ck_mem.ck_znQ[0] == null) {\r
+             N_VDestroy(ck_mem.ck_zn[0]);\r
+             N_VDestroy(ck_mem.ck_zn[1]);\r
+             ck_mem = null;\r
+             return(null);\r
+           }\r
 \r
-                 /* Load derivative */\r
+           N_VScale_Serial(ONE, cv_mem.cv_znQ[0], ck_mem.ck_znQ[0]);\r
 \r
-                 if (cv_mem.cv_nst == 0) {\r
+         }\r
 \r
-                   /**  retval = */ \r
-                         if (testMode) {\r
-                                 fTestMode(cv_mem.cv_tn, content.y, content.yd, cv_mem.cv_user_data);\r
-                         }\r
-                         else {\r
-                             f(cv_mem.cv_tn, content.y, content.yd, cv_mem.cv_user_data);\r
-                         }\r
+         /* Do we need to carry sensitivities? */\r
+         ck_mem.ck_sensi = cv_mem.cv_sensi;\r
 \r
-                   if (ca_mem.ca_IMstoreSensi) {\r
-                     /* retval = */ cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, content.y, content.yd,\r
-                                               content.yS, content.ySd,\r
-                                               cv_mem.cv_tempv, cv_mem.cv_ftemp);\r
-                   }\r
+         if (ck_mem.ck_sensi) {\r
 \r
-                 } else {\r
+           ck_mem.ck_Ns = cv_mem.cv_Ns;\r
 \r
-                   N_VScale_Serial(ONE/cv_mem.cv_h, cv_mem.cv_zn[1], content.yd);\r
+           ck_mem.ck_znS[0] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+           if (ck_mem.ck_znS[0] == null) {\r
+             N_VDestroy(ck_mem.ck_zn[0]);\r
+             N_VDestroy(ck_mem.ck_zn[1]);\r
+             if (ck_mem.ck_quadr) N_VDestroy(ck_mem.ck_znQ[0]);\r
+             ck_mem = null;\r
+             return(null);\r
+           }\r
 \r
-                   if (ca_mem.ca_IMstoreSensi) {\r
-                     for (is=0; is<cv_mem.cv_Ns; is++) \r
-                       N_VScale_Serial(ONE/cv_mem.cv_h, cv_mem.cv_znS[1][is], content.ySd[is]);\r
-                   }\r
+           for (is=0; is<cv_mem.cv_Ns; is++)\r
+             N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], ck_mem.ck_znS[0][is]);\r
 \r
-                 }\r
+         }\r
 \r
-                 return(0);\r
-               }\r
+         /* Do we need to carry quadrature sensitivities? */\r
+         ck_mem.ck_quadr_sensi = cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS;\r
+\r
+         if (ck_mem.ck_quadr_sensi) {\r
+           ck_mem.ck_znQS[0] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
+           if (ck_mem.ck_znQS[0] == null) {\r
+             N_VDestroy(ck_mem.ck_zn[0]);\r
+             N_VDestroy(ck_mem.ck_zn[1]);\r
+             if (ck_mem.ck_quadr) N_VDestroy(ck_mem.ck_znQ[0]);\r
+             N_VDestroyVectorArray_Serial(ck_mem.ck_znS[0], cv_mem.cv_Ns);\r
+             ck_mem = null;\r
+             return(null);\r
+           }\r
+           \r
+           for (is=0; is<cv_mem.cv_Ns; is++)\r
+             N_VScale_Serial(ONE, cv_mem.cv_znQS[0][is], ck_mem.ck_znQS[0][is]);\r
 \r
-               /*\r
-                * CVApolynomialStorePnt ( -> IMstore )\r
-                *\r
-                * This routine stores a new point y in the structure d for use\r
-                * in the Polynomial interpolation.\r
-                * Note that the time is already stored.\r
-                */\r
+         }\r
 \r
-               private int CVApolynomialStorePnt(CVodeMemRec cv_mem, DtpntMemRec d)\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 PolynomialDataMemRec content;\r
-                 int is;\r
+         /* Next in list */\r
+         ck_mem.ck_next  = null;\r
 \r
-                 ca_mem = cv_mem.cv_adj_mem;\r
+         return(ck_mem);\r
+       }\r
+       \r
+       /*\r
+        * CVAhermiteMalloc\r
+        *\r
+        * This routine allocates memory for storing information at all\r
+        * intermediate points between two consecutive check points. \r
+        * This data is then used to interpolate the forward solution \r
+        * at any other time.\r
+        */\r
 \r
-                 content = d.polynomialContent;\r
+       private boolean CVAhermiteMalloc(CVodeMemRec cv_mem)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         DtpntMemRec dt_mem[];\r
+         HermiteDataMemRec content;\r
+         int i, ii=0;\r
+         boolean allocOK;\r
 \r
-                 N_VScale_Serial(ONE, cv_mem.cv_zn[0], content.y);\r
+         allocOK = true;\r
 \r
-                 if (ca_mem.ca_IMstoreSensi) {\r
-                   for (is=0; is<cv_mem.cv_Ns; is++) \r
-                     N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], content.yS[is]);\r
-                 }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-                 content.order = cv_mem.cv_qu;\r
+         /* Allocate space for the vectors ytmp and yStmp */\r
 \r
-                 return(0);\r
-               }\r
-               \r
-               /*\r
-                * CVAckpntNew\r
-                *\r
-                * This routine allocates space for a new check point and sets \r
-                * its data from current values in cv_mem.\r
-                */\r
-\r
-               private CkpntMemRec CVAckpntNew(CVodeMemRec cv_mem)\r
-               {\r
-                 CkpntMemRec ck_mem;\r
-                 int j, jj, is, qmax; \r
-\r
-                 /* Allocate space for ckdata */\r
-                 ck_mem = null;\r
-                 ck_mem = new CkpntMemRec();\r
-\r
-                 /* Set cv_next to NULL */\r
-                 ck_mem.ck_next = null;\r
-\r
-                 /* Test if we need to allocate space for the last zn.\r
-                  * NOTE: zn(qmax) may be needed for a hot restart, if an order\r
-                  * increase is deemed necessary at the first step after a check point */\r
-                 qmax = cv_mem.cv_qmax;\r
-                 ck_mem.ck_zqm = (cv_mem.cv_q < qmax) ? qmax : 0;\r
-\r
-                 for (j=0; j<=cv_mem.cv_q; j++) {\r
-                   ck_mem.ck_zn[j] = N_VClone(cv_mem.cv_tempv);\r
-                   if (ck_mem.ck_zn[j] == null) {\r
-                     for (jj=0; jj<j; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                     ck_mem = null;\r
-                     return(null);\r
-                   }\r
-                 }\r
+         ca_mem.ca_ytmp = N_VClone(cv_mem.cv_tempv);\r
+         if (ca_mem.ca_ytmp == null) {\r
+           return(false);\r
+         }\r
 \r
-                 if (cv_mem.cv_q < qmax) {\r
-                   ck_mem.ck_zn[qmax] = N_VClone(cv_mem.cv_tempv);\r
-                   if (ck_mem.ck_zn[qmax] == null) {\r
-                     for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                     ck_mem = null;\r
-                     return(null);\r
-                   }\r
-                 }\r
+         if (ca_mem.ca_IMstoreSensi) {\r
+           ca_mem.ca_yStmp = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+           if (ca_mem.ca_yStmp == null) {\r
+             N_VDestroy(ca_mem.ca_ytmp);\r
+             return(false);\r
+           }\r
+         }\r
 \r
-                 /* Test if we need to carry quadratures */\r
-                 ck_mem.ck_quadr = cv_mem.cv_quadr && cv_mem.cv_errconQ;\r
-\r
-                 if (ck_mem.ck_quadr) {\r
-\r
-                   for (j=0; j<=cv_mem.cv_q; j++) {\r
-                     ck_mem.ck_znQ[j] = N_VClone(cv_mem.cv_tempvQ);\r
-                     if(ck_mem.ck_znQ[j] == null) {\r
-                       for (jj=0; jj<j; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
-                       if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; j++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
-\r
-                   if (cv_mem.cv_q < qmax) {\r
-                     ck_mem.ck_znQ[qmax] = N_VClone(cv_mem.cv_tempvQ);\r
-                     if (ck_mem.ck_znQ[qmax] == null) {\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
-                       N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
+         /* Allocate space for the content field of the dt structures */\r
 \r
-                 }\r
+         dt_mem = ca_mem.dt_mem;\r
 \r
-                 /* Test if we need to carry sensitivities */\r
-                 ck_mem.ck_sensi = cv_mem.cv_sensi;\r
-\r
-                 if (ck_mem.ck_sensi) {\r
-\r
-                   ck_mem.ck_Ns = cv_mem.cv_Ns;\r
-\r
-                   for (j=0; j<=cv_mem.cv_q; j++) {\r
-                     ck_mem.ck_znS[j] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                     if (ck_mem.ck_znS[j] == null) {\r
-                       for (jj=0; jj<j; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
-                       if (ck_mem.ck_quadr) {\r
-                         if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_znQ[qmax]);\r
-                         for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
-                       }\r
-                       if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
-\r
-                   if ( cv_mem.cv_q < qmax) {\r
-                     ck_mem.ck_znS[qmax] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
-                     if (ck_mem.ck_znS[qmax] == null) {\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
-                       if (ck_mem.ck_quadr) {\r
-                         N_VDestroy(ck_mem.ck_znQ[qmax]);\r
-                         for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
-                       }\r
-                       N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
+         for (i=0; i<=ca_mem.ca_nsteps; i++) {\r
 \r
-                 }\r
+           content = null;\r
+           content = new HermiteDataMemRec();\r
 \r
-                 /* Test if we need to carry quadrature sensitivities */\r
-                 ck_mem.ck_quadr_sensi = cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS;\r
-\r
-                 if (ck_mem.ck_quadr_sensi) {\r
-\r
-                   for (j=0; j<=cv_mem.cv_q; j++) {\r
-                     ck_mem.ck_znQS[j] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
-                     if (ck_mem.ck_znQS[j] == null) {\r
-                       for (jj=0; jj<j; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znQS[jj], cv_mem.cv_Ns);\r
-                       if (cv_mem.cv_q < qmax) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[qmax], cv_mem.cv_Ns);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
-                       if (ck_mem.ck_quadr) {\r
-                         if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_znQ[qmax]);\r
-                         for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
-                       }\r
-                       if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
-\r
-                   if ( cv_mem.cv_q < qmax) {\r
-                     ck_mem.ck_znQS[qmax] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
-                     if (ck_mem.ck_znQS[qmax] == null) {\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znQS[jj], cv_mem.cv_Ns);\r
-                       N_VDestroyVectorArray_Serial(ck_mem.ck_znS[qmax], cv_mem.cv_Ns);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
-                       if (ck_mem.ck_quadr) {\r
-                         N_VDestroy(ck_mem.ck_znQ[qmax]);\r
-                         for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       }\r
-                       N_VDestroy(ck_mem.ck_zn[qmax]);\r
-                       for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
-                       ck_mem = null;\r
-                       return(null);\r
-                     }\r
-                   }\r
+           content.y = N_VClone(cv_mem.cv_tempv);\r
+           if (content.y == null) {\r
+             content = null;\r
+             ii = i;\r
+             allocOK = false;\r
+             break;\r
+           }\r
 \r
-                 }\r
+           content.yd = N_VClone(cv_mem.cv_tempv);\r
+           if (content.yd == null) {\r
+             N_VDestroy(content.y);\r
+             content = null;\r
+             ii = i;\r
+             allocOK = false;\r
+             break;\r
+           }\r
 \r
-                 /* Load check point data from cv_mem */\r
+           if (ca_mem.ca_IMstoreSensi) {\r
 \r
-                 for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_zn[j], ck_mem.ck_zn[j]);\r
-                 if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_zn[qmax], ck_mem.ck_zn[qmax]);\r
+             content.yS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+             if (content.yS == null) {\r
+               N_VDestroy(content.y);\r
+               N_VDestroy(content.yd);\r
+               content = null;\r
+               ii = i;\r
+               allocOK = false;\r
+               break;\r
+             }\r
 \r
-                 if (ck_mem.ck_quadr) {\r
-                   for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znQ[j], ck_mem.ck_znQ[j]);\r
-                   if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znQ[qmax], ck_mem.ck_znQ[qmax]);\r
-                 }\r
+             content.ySd = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+             if (content.ySd == null) {\r
+               N_VDestroy(content.y);\r
+               N_VDestroy(content.yd);\r
+               N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
+               content = null;\r
+               ii = i;\r
+               allocOK = false;\r
+               break;\r
+             }\r
+             \r
+           }\r
+           \r
+           dt_mem[i].content = CV_HERMITE;\r
+           dt_mem[i].hermiteContent = content;\r
 \r
-                 if (ck_mem.ck_sensi) {\r
-                   for (is=0; is<cv_mem.cv_Ns; is++) {\r
-                     for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znS[j][is], ck_mem.ck_znS[j][is]);\r
-                     if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znS[qmax][is], ck_mem.ck_znS[qmax][is]);\r
-                   }\r
-                 }\r
+         } \r
 \r
-                 if (ck_mem.ck_quadr_sensi) {\r
-                   for (is=0; is<cv_mem.cv_Ns; is++) {\r
-                     for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znQS[j][is], ck_mem.ck_znQS[j][is]);\r
-                     if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znQS[qmax][is], ck_mem.ck_znQS[qmax][is]);\r
-                   }\r
-                 }\r
+         /* If an error occurred, deallocate and return */\r
 \r
-                 for (j=0; j<=L_MAX; j++)        ck_mem.ck_tau[j] = cv_mem.cv_tau[j];\r
-                 for (j=0; j<=NUM_TESTS; j++)    ck_mem.ck_tq[j] = cv_mem.cv_tq[j];\r
-                 for (j=0; j<=cv_mem.cv_q; j++) ck_mem.ck_l[j] = cv_mem.cv_l[j];\r
-                 ck_mem.ck_nst       = cv_mem.cv_nst;\r
-                 ck_mem.ck_tretlast  = cv_mem.cv_tretlast;\r
-                 ck_mem.ck_q         = cv_mem.cv_q;\r
-                 ck_mem.ck_qprime    = cv_mem.cv_qprime;\r
-                 ck_mem.ck_qwait     = cv_mem.cv_qwait;\r
-                 ck_mem.ck_L         = cv_mem.cv_L;\r
-                 ck_mem.ck_gammap    = cv_mem.cv_gammap;\r
-                 ck_mem.ck_h         = cv_mem.cv_h;\r
-                 ck_mem.ck_hprime    = cv_mem.cv_hprime;\r
-                 ck_mem.ck_hscale    = cv_mem.cv_hscale;\r
-                 ck_mem.ck_eta       = cv_mem.cv_eta;\r
-                 ck_mem.ck_etamax    = cv_mem.cv_etamax;\r
-                 ck_mem.ck_t0        = cv_mem.cv_tn;\r
-                 ck_mem.ck_saved_tq5 = cv_mem.cv_saved_tq5;\r
-\r
-                 return(ck_mem);\r
-               }\r
-               \r
-               /* \r
-                * CVodeGetQuad\r
-                *\r
-                * This routine extracts quadrature solution into yQout at the\r
-                * time which CVode returned the solution.\r
-                * This is just a wrapper that calls CVodeGetQuadDky with k=0.\r
-                */\r
-                \r
-               private int CVodeGetQuad(CVodeMemRec cv_mem, double tret[], NVector yQout)\r
-               {\r
-                 int flag;\r
-\r
-                 if (cv_mem == null) {\r
-                   cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuad", MSGCV_NO_MEM);\r
-                   return(CV_MEM_NULL);\r
-                 } \r
-\r
-                 tret[0] = cv_mem.cv_tretlast;\r
-                 \r
-                 flag = CVodeGetQuadDky(cv_mem,cv_mem.cv_tretlast,0,yQout);\r
-\r
-                 return(flag);\r
-               }\r
-               \r
-               /*\r
-                * CVodeGetQuadDky\r
-                *\r
-                * CVodeQuadDky computes the kth derivative of the yQ function at\r
-                * time t, where tn-hu <= t <= tn, tn denotes the current         \r
-                * internal time reached, and hu is the last internal step size   \r
-                * successfully used by the solver. The user may request \r
-                * k=0, 1, ..., qu, where qu is the current order. \r
-                * The derivative vector is returned in dky. This vector \r
-                * must be allocated by the caller. It is only legal to call this         \r
-                * function after a successful return from CVode with quadrature\r
-                * computation enabled.\r
-                */\r
-\r
-               private int CVodeGetQuadDky(CVodeMemRec cv_mem, double t, int k, NVector dkyQ)\r
-               { \r
-                 double s, c, r;\r
-                 double tfuzz, tp, tn1;\r
-                 int i, j;\r
-                 \r
-                 /* Check all inputs for legality */\r
-                 \r
-                 if (cv_mem == null) {\r
-                   cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuadDky", MSGCV_NO_MEM);\r
-                   return(CV_MEM_NULL);\r
-                 } \r
-\r
-                 if(cv_mem.cv_quadr != true) {\r
-                   cvProcessError(cv_mem, CV_NO_QUAD, "CVODES", "CVodeGetQuadDky", MSGCV_NO_QUAD);\r
-                   return(CV_NO_QUAD);\r
-                 }\r
+         if (!allocOK) {\r
 \r
-                 if (dkyQ == null) {\r
-                   cvProcessError(cv_mem, CV_BAD_DKY, "CVODES", "CVodeGetQuadDky", MSGCV_NULL_DKY);\r
-                   return(CV_BAD_DKY);\r
-                 }\r
-                 \r
-                 if ((k < 0) || (k > cv_mem.cv_q)) {\r
-                   cvProcessError(cv_mem, CV_BAD_K, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_K);\r
-                   return(CV_BAD_K);\r
-                 }\r
-                 \r
-                 /* Allow for some slack */\r
-                 tfuzz = FUZZ_FACTOR * cv_mem.cv_uround *\r
-                   (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_hu));\r
-                 if (cv_mem.cv_hu < ZERO) tfuzz = -tfuzz;\r
-                 tp = cv_mem.cv_tn - cv_mem.cv_hu - tfuzz;\r
-                 tn1 = cv_mem.cv_tn + tfuzz;\r
-                 if ((t-tp)*(t-tn1) > ZERO) {\r
-                   cvProcessError(cv_mem, CV_BAD_T, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_T);\r
-                   return(CV_BAD_T);\r
-                 }\r
-                 \r
-                 /* Sum the differentiated interpolating polynomial */\r
-                 \r
-                 s = (t - cv_mem.cv_tn) / cv_mem.cv_h;\r
-                 for (j=cv_mem.cv_q; j >= k; j--) {\r
-                   c = ONE;\r
-                   for (i=j; i >= j-k+1; i--) c *= i;\r
-                   if (j == cv_mem.cv_q) {\r
-                     N_VScale_Serial(c, cv_mem.cv_znQ[cv_mem.cv_q], dkyQ);\r
-                   } else {\r
-                     N_VLinearSum_Serial(c, cv_mem.cv_znQ[j], s, dkyQ, dkyQ);\r
-                   }\r
-                 }\r
-                 if (k == 0) return(CV_SUCCESS);\r
-                 r = Math.pow(cv_mem.cv_h, -k);\r
-                 N_VScale_Serial(r, dkyQ, dkyQ);\r
-                 return(CV_SUCCESS);\r
-                 \r
-               }\r
+           N_VDestroy(ca_mem.ca_ytmp);\r
 \r
-               class CVadjCheckPointRec {\r
-                         CkpntMemRec my_addr;\r
-                         CkpntMemRec next_addr;\r
-                         double t0;\r
-                         double t1;\r
-                         long nstep;\r
-                         int order;\r
-                         double step;\r
-                       };\r
-\r
-                       /*\r
-                        * CVodeGetAdjCheckPointsInfo\r
-                        *\r
-                        * This routine loads an array of nckpnts structures of type CVadjCheckPointRec.\r
-                        * The user must allocate space for ckpnt.\r
-                        */\r
-\r
-                       private int CVodeGetAdjCheckPointsInfo(CVodeMemRec cv_mem, CVadjCheckPointRec ckpnt[])\r
-                       {\r
-                         CVadjMemRec ca_mem;\r
-                         CkpntMemRec ck_mem;\r
-                         int i;\r
-\r
-                         /* Check if cvode_mem exists */\r
-                         if (cv_mem == null) {\r
-                           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_MEM);\r
-                           return(CV_MEM_NULL);\r
-                         }\r
+           if (ca_mem.ca_IMstoreSensi) {\r
+             N_VDestroyVectorArray_Serial(ca_mem.ca_yStmp, cv_mem.cv_Ns);\r
+           }\r
 \r
-                         /* Was ASA initialized? */\r
-                         if (cv_mem.cv_adjMallocDone == false) {\r
-                           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_ADJ);\r
-                           return(CV_NO_ADJ);\r
-                         } \r
-                         ca_mem = cv_mem.cv_adj_mem;\r
+           for (i=0; i<ii; i++) {\r
+             content = (dt_mem[i].hermiteContent);\r
+             N_VDestroy(content.y);\r
+             N_VDestroy(content.yd);\r
+             if (ca_mem.ca_IMstoreSensi) {\r
+               N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
+               N_VDestroyVectorArray_Serial(content.ySd, cv_mem.cv_Ns);\r
+             }\r
+             dt_mem[i].hermiteContent = null;\r
+           }\r
 \r
-                         ck_mem = ca_mem.ck_mem;\r
+         }\r
 \r
-                         i = 0;\r
+         return(allocOK);\r
+       }\r
 \r
-                         while (ck_mem != null) {\r
+       \r
+       /*\r
+        * CVApolynomialMalloc\r
+        *\r
+        * This routine allocates memory for storing information at all\r
+        * intermediate points between two consecutive check points. \r
+        * This data is then used to interpolate the forward solution \r
+        * at any other time.\r
+        */\r
 \r
-                           ckpnt[i].my_addr = ck_mem;\r
-                           ckpnt[i].next_addr = ck_mem.ck_next;\r
-                           ckpnt[i].t0 = ck_mem.ck_t0;\r
-                           ckpnt[i].t1 = ck_mem.ck_t1;\r
-                           ckpnt[i].nstep = ck_mem.ck_nst;\r
-                           ckpnt[i].order = ck_mem.ck_q;\r
-                           ckpnt[i].step = ck_mem.ck_h;\r
+       private boolean CVApolynomialMalloc(CVodeMemRec cv_mem)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         DtpntMemRec dt_mem[];\r
+         PolynomialDataMemRec content;\r
+         int i, ii=0;\r
+         boolean allocOK;\r
 \r
-                           ck_mem = ck_mem.ck_next;\r
-                           i++;\r
+         allocOK = true;\r
 \r
-                         }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-                         return(CV_SUCCESS);\r
+         /* Allocate space for the vectors ytmp and yStmp */\r
 \r
-                       }\r
-                       \r
-                       private int CVodeCreateB(CVodeMemRec cv_mem, int lmmB, int iterB, int which[])\r
-                       {\r
-                         CVadjMemRec ca_mem;\r
-                         CVodeBMemRec new_cvB_mem;\r
-                         CVodeMemRec cvodeB_mem;\r
+         ca_mem.ca_ytmp = N_VClone(cv_mem.cv_tempv);\r
+         if (ca_mem.ca_ytmp == null) {\r
+           return(false);\r
+         }\r
 \r
-                         /* Check if cvode_mem exists */\r
-                         if (cv_mem == null) {\r
-                           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeCreateB", MSGCV_NO_MEM);\r
-                           return(CV_MEM_NULL);\r
-                         }\r
+         if (ca_mem.ca_IMstoreSensi) {\r
+           ca_mem.ca_yStmp = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+           if (ca_mem.ca_yStmp == null) {\r
+             N_VDestroy(ca_mem.ca_ytmp);\r
+             return(false);\r
+           }\r
+         }\r
 \r
-                         /* Was ASA initialized? */\r
-                         if (cv_mem.cv_adjMallocDone == false) {\r
-                           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeCreateB", MSGCV_NO_ADJ);\r
-                           return(CV_NO_ADJ);\r
-                         }\r
-                         ca_mem = cv_mem.cv_adj_mem;\r
+         /* Allocate space for the content field of the dt structures */\r
 \r
-                         /* Allocate space for new CVodeBMem object */\r
+         dt_mem = ca_mem.dt_mem;\r
 \r
-                         new_cvB_mem = null;\r
-                         new_cvB_mem = new CVodeBMemRec();\r
+         for (i=0; i<=ca_mem.ca_nsteps; i++) {\r
 \r
-                         /* Create and set a new CVODES object for the backward problem */\r
+           content = null;\r
+           content = new PolynomialDataMemRec();\r
 \r
-                         cvodeB_mem = CVodeCreate(lmmB, iterB);\r
-                         if (cvodeB_mem == null) {\r
-                           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeCreateB", MSGCV_MEM_FAIL);\r
-                           return(CV_MEM_FAIL);\r
-                         }\r
+           content.y = N_VClone(cv_mem.cv_tempv);\r
+           if (content.y == null) {\r
+             content = null;\r
+             ii = i;\r
+             allocOK = false;\r
+             break;\r
+           }\r
 \r
-                         UserData userData = new UserData();\r
-                         userData.memRec = cv_mem;\r
-                         cvodeB_mem.cv_user_data = userData;\r
+           if (ca_mem.ca_IMstoreSensi) {\r
 \r
-                         cvodeB_mem.cv_mxhnil = -1;\r
+             content.yS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+             if (content.yS == null) {\r
+               N_VDestroy(content.y);\r
+               content = null;\r
+               ii = i;\r
+               allocOK = false;\r
+               break;\r
+             }\r
 \r
-                         //CVodeSetErrHandlerFn(cvodeB_mem, cv_mem->cv_ehfun, cv_mem->cv_eh_data);\r
-                         //cvodeB_mem.cv_ehfun = cv_mem.cv_ehfun;\r
-                         cvodeB_mem.cv_eh_data = cv_mem.cv_eh_data;\r
+           }\r
 \r
-                         //CVodeSetErrFile(cvodeB_mem, cv_mem->cv_errfp);\r
-                         cvodeB_mem.cv_errfp = cv_mem.cv_errfp;\r
+           dt_mem[i].content = CV_POLYNOMIAL;\r
+           dt_mem[i].polynomialContent = content;\r
 \r
-                         /* Set/initialize fields in the new CVodeBMem object, new_cvB_mem */\r
+         } \r
 \r
-                         new_cvB_mem.cv_index   = ca_mem.ca_nbckpbs;\r
+         /* If an error occurred, deallocate and return */\r
 \r
-                         new_cvB_mem.cv_mem     = cvodeB_mem;\r
+         if (!allocOK) {\r
 \r
-                         //new_cvB_mem.cv_f       = null;\r
-                         //new_cvB_mem.cv_fs      = null;\r
+           N_VDestroy(ca_mem.ca_ytmp);\r
 \r
-                         //new_cvB_mem.cv_fQ      = null;\r
-                         //new_cvB_mem.cv_fQs     = null;\r
+           if (ca_mem.ca_IMstoreSensi) {\r
+             N_VDestroyVectorArray_Serial(ca_mem.ca_yStmp, cv_mem.cv_Ns);\r
+           }\r
 \r
-                         new_cvB_mem.cv_user_data  = null;\r
+           for (i=0; i<ii; i++) {\r
+             content = (dt_mem[i].polynomialContent);\r
+             N_VDestroy(content.y);\r
+             if (ca_mem.ca_IMstoreSensi) {\r
+               N_VDestroyVectorArray_Serial(content.yS, cv_mem.cv_Ns);\r
+             }\r
+             dt_mem[i].polynomialContent = null;\r
+           }\r
 \r
-                         //new_cvB_mem.cv_lmem    = null;\r
-                         //new_cvB_mem.cv_lfree   = null;\r
-                         //new_cvB_mem.cv_pmem    = null;\r
-                         //new_cvB_mem.cv_pfree   = null;\r
+         }\r
 \r
-                         new_cvB_mem.cv_y       = null;\r
+         return(allocOK);\r
 \r
-                         new_cvB_mem.cv_f_withSensi = false;\r
-                         new_cvB_mem.cv_fQ_withSensi = false;\r
+       }\r
+       \r
+       /*\r
+        * CVAhermiteStorePnt ( -> IMstore )\r
+        *\r
+        * This routine stores a new point (y,yd) in the structure d for use\r
+        * in the cubic Hermite interpolation.\r
+        * Note that the time is already stored.\r
+        */\r
 \r
-                         /* Attach the new object to the linked list cvB_mem */\r
+       private int CVAhermiteStorePnt(CVodeMemRec cv_mem, DtpntMemRec d)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         HermiteDataMemRec content;\r
+         int is;\r
+         /* int retval; */\r
 \r
-                         new_cvB_mem.cv_next = ca_mem.cvB_mem;\r
-                         ca_mem.cvB_mem = new_cvB_mem;\r
-                         \r
-                         /* Return the index of the newly created CVodeBMem object.\r
-                          * This must be passed to CVodeInitB and to other ***B \r
-                          * functions to set optional inputs for this backward problem */\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-                         which[0] = ca_mem.ca_nbckpbs;\r
+         content = d.hermiteContent;\r
 \r
-                         ca_mem.ca_nbckpbs++;\r
+         /* Load solution */\r
 \r
-                         return(CV_SUCCESS);\r
-                       }\r
+         N_VScale_Serial(ONE, cv_mem.cv_zn[0], content.y);\r
+         \r
+         if (ca_mem.ca_IMstoreSensi) {\r
+           for (is=0; is<cv_mem.cv_Ns; is++) \r
+             N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], content.yS[is]);\r
+         }\r
 \r
-                       private int CVodeInitB(CVodeMemRec cv_mem, int which, \r
-                              int fB,\r
-                              double tB0, NVector yB0)\r
-               {\r
-                 CVadjMemRec ca_mem;\r
-                 CVodeBMemRec cvB_mem;\r
-                 CVodeMemRec cvodeB_mem;\r
-                 int flag;\r
+         /* Load derivative */\r
 \r
-                 /* Check if cvode_mem exists */\r
+         if (cv_mem.cv_nst == 0) {\r
 \r
-                 if (cv_mem == null) {\r
-                   cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeInitB", MSGCV_NO_MEM);\r
-                   return(CV_MEM_NULL);\r
+           /**  retval = */ \r
+                 if (testMode) {\r
+                         fTestMode(cv_mem.cv_tn, content.y, content.yd, cv_mem.cv_user_data);\r
                  }\r
-\r
-                 /* Was ASA initialized? */\r
-\r
-                 if (cv_mem.cv_adjMallocDone == false) {\r
-                   cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeInitB", MSGCV_NO_ADJ);\r
-                   return(CV_NO_ADJ);\r
-                 } \r
-                 ca_mem = cv_mem.cv_adj_mem;\r
-\r
-                 /* Check the value of which */\r
-\r
-                 if ( which >= ca_mem.ca_nbckpbs ) {\r
-                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeInitB", MSGCV_BAD_WHICH);\r
-                   return(CV_ILL_INPUT);\r
+                 else {\r
+                     f(cv_mem.cv_tn, content.y, content.yd, cv_mem.cv_user_data);\r
                  }\r
 \r
-                 /* Find the CVodeBMem entry in the linked list corresponding to which */\r
-\r
-                 cvB_mem = ca_mem.cvB_mem;\r
-                 while (cvB_mem != null) {\r
-                   if ( which == cvB_mem.cv_index ) break;\r
-                   cvB_mem = cvB_mem.cv_next;\r
-                 }\r
+           if (ca_mem.ca_IMstoreSensi) {\r
+             /* retval = */ cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, content.y, content.yd,\r
+                                       content.yS, content.ySd,\r
+                                       cv_mem.cv_tempv, cv_mem.cv_ftemp);\r
+           }\r
 \r
-                 cvodeB_mem = cvB_mem.cv_mem;\r
-                 \r
-                 /* Allocate and set the CVODES object */\r
+         } else {\r
 \r
-                 flag = CVodeInit(cvodeB_mem, CVArhs_select, tB0, yB0);\r
+           N_VScale_Serial(ONE/cv_mem.cv_h, cv_mem.cv_zn[1], content.yd);\r
 \r
-                 if (flag != CV_SUCCESS) return(flag);\r
+           if (ca_mem.ca_IMstoreSensi) {\r
+             for (is=0; is<cv_mem.cv_Ns; is++) \r
+               N_VScale_Serial(ONE/cv_mem.cv_h, cv_mem.cv_znS[1][is], content.ySd[is]);\r
+           }\r
 \r
-                 /* Copy fB function in cvB_mem */\r
+         }\r
 \r
-                 cvB_mem.cv_f_withSensi = false;\r
-                 //cvB_mem.cv_f = fB;\r
+         return(0);\r
+       }\r
 \r
-                 /* Allocate space and initialize the y Nvector in cvB_mem */\r
+       /*\r
+        * CVApolynomialStorePnt ( -> IMstore )\r
+        *\r
+        * This routine stores a new point y in the structure d for use\r
+        * in the Polynomial interpolation.\r
+        * Note that the time is already stored.\r
+        */\r
 \r
-                 cvB_mem.cv_t0 = tB0;\r
-                 cvB_mem.cv_y = N_VClone(yB0);\r
-                 N_VScale_Serial(ONE, yB0, cvB_mem.cv_y);\r
+       private int CVApolynomialStorePnt(CVodeMemRec cv_mem, DtpntMemRec d)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         PolynomialDataMemRec content;\r
+         int is;\r
 \r
-                 return(CV_SUCCESS);\r
-               }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-                       /*\r
-                        * CVArhs\r
-                        *\r
-                        * This routine interfaces to the CVRhsFnB (or CVRhsFnBS) routine \r
-                        * provided by the user.\r
-                        */\r
+         content = d.polynomialContent;\r
 \r
+         N_VScale_Serial(ONE, cv_mem.cv_zn[0], content.y);\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
+         if (ca_mem.ca_IMstoreSensi) {\r
+           for (is=0; is<cv_mem.cv_Ns; is++) \r
+             N_VScale_Serial(ONE, cv_mem.cv_znS[0][is], content.yS[is]);\r
+         }\r
 \r
-                         ca_mem = cv_mem.cv_adj_mem;\r
+         content.order = cv_mem.cv_qu;\r
 \r
-                         cvB_mem = ca_mem.ca_bckpbCrt;\r
+         return(0);\r
+       }\r
+       \r
+       /*\r
+        * CVAckpntNew\r
+        *\r
+        * This routine allocates space for a new check point and sets \r
+        * its data from current values in cv_mem.\r
+        */\r
 \r
-                         /* Get forward solution from interpolation */\r
+       private CkpntMemRec CVAckpntNew(CVodeMemRec cv_mem)\r
+       {\r
+         CkpntMemRec ck_mem;\r
+         int j, jj, is, qmax; \r
+\r
+         /* Allocate space for ckdata */\r
+         ck_mem = null;\r
+         ck_mem = new CkpntMemRec();\r
+\r
+         /* Set cv_next to NULL */\r
+         ck_mem.ck_next = null;\r
+\r
+         /* Test if we need to allocate space for the last zn.\r
+          * NOTE: zn(qmax) may be needed for a hot restart, if an order\r
+          * increase is deemed necessary at the first step after a check point */\r
+         qmax = cv_mem.cv_qmax;\r
+         ck_mem.ck_zqm = (cv_mem.cv_q < qmax) ? qmax : 0;\r
+\r
+         for (j=0; j<=cv_mem.cv_q; j++) {\r
+           ck_mem.ck_zn[j] = N_VClone(cv_mem.cv_tempv);\r
+           if (ck_mem.ck_zn[j] == null) {\r
+             for (jj=0; jj<j; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+             ck_mem = null;\r
+             return(null);\r
+           }\r
+         }\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
-                                 if (ca_mem.ca_IMget == CVAhermiteGetY_select) {\r
-                                 flag = CVAhermiteGetY(cv_mem, t, ca_mem.ca_ytmp, null);\r
-                                 }\r
+         if (cv_mem.cv_q < qmax) {\r
+           ck_mem.ck_zn[qmax] = N_VClone(cv_mem.cv_tempv);\r
+           if (ck_mem.ck_zn[qmax] == null) {\r
+             for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+             ck_mem = null;\r
+             return(null);\r
+           }\r
+         }\r
 \r
-                         if (flag != CV_SUCCESS) {\r
-                           cvProcessError(cv_mem, -1, "CVODEA", "CVArhs", MSGCV_BAD_TINTERP, t);\r
-                           return(-1);\r
-                         }\r
+         /* Test if we need to carry quadratures */\r
+         ck_mem.ck_quadr = cv_mem.cv_quadr && cv_mem.cv_errconQ;\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
-                        * 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
-                       private 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
+         if (ck_mem.ck_quadr) {\r
 \r
-                         /* Extract stuff from the appropriate data points */\r
+           for (j=0; j<=cv_mem.cv_q; j++) {\r
+             ck_mem.ck_znQ[j] = N_VClone(cv_mem.cv_tempvQ);\r
+             if(ck_mem.ck_znQ[j] == null) {\r
+               for (jj=0; jj<j; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
+               if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; j++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\r
+           }\r
 \r
-                         t0 = dt_mem[(int)(indx[0]-1)].t;\r
-                         t1 = dt_mem[(int)indx[0]].t;\r
-                         delta = t1 - t0;\r
+           if (cv_mem.cv_q < qmax) {\r
+             ck_mem.ck_znQ[qmax] = N_VClone(cv_mem.cv_tempvQ);\r
+             if (ck_mem.ck_znQ[qmax] == null) {\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
+               N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\r
+           }\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
 \r
-                         if (newpoint[0]) {\r
-                           \r
-                           /* Recompute Y0 and Y1 */\r
+         /* Test if we need to carry sensitivities */\r
+         ck_mem.ck_sensi = cv_mem.cv_sensi;\r
+\r
+         if (ck_mem.ck_sensi) {\r
+\r
+           ck_mem.ck_Ns = cv_mem.cv_Ns;\r
+\r
+           for (j=0; j<=cv_mem.cv_q; j++) {\r
+             ck_mem.ck_znS[j] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+             if (ck_mem.ck_znS[j] == null) {\r
+               for (jj=0; jj<j; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
+               if (ck_mem.ck_quadr) {\r
+                 if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_znQ[qmax]);\r
+                 for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
+               }\r
+               if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\r
+           }\r
 \r
-                           content1 = dt_mem[(int)indx[0]].hermiteContent;\r
+           if ( cv_mem.cv_q < qmax) {\r
+             ck_mem.ck_znS[qmax] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempv);\r
+             if (ck_mem.ck_znS[qmax] == null) {\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
+               if (ck_mem.ck_quadr) {\r
+                 N_VDestroy(ck_mem.ck_znQ[qmax]);\r
+                 for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
+               }\r
+               N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\r
+           }\r
 \r
-                           y1  = content1.y;\r
-                           yd1 = content1.yd;\r
+         }\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
+         /* Test if we need to carry quadrature sensitivities */\r
+         ck_mem.ck_quadr_sensi = cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS;\r
+\r
+         if (ck_mem.ck_quadr_sensi) {\r
+\r
+           for (j=0; j<=cv_mem.cv_q; j++) {\r
+             ck_mem.ck_znQS[j] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
+             if (ck_mem.ck_znQS[j] == null) {\r
+               for (jj=0; jj<j; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znQS[jj], cv_mem.cv_Ns);\r
+               if (cv_mem.cv_q < qmax) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[qmax], cv_mem.cv_Ns);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
+               if (ck_mem.ck_quadr) {\r
+                 if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_znQ[qmax]);\r
+                 for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_znQ[jj]);\r
+               }\r
+               if (cv_mem.cv_q < qmax) N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\r
+           }\r
 \r
+           if ( cv_mem.cv_q < qmax) {\r
+             ck_mem.ck_znQS[qmax] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, cv_mem.cv_tempvQ);\r
+             if (ck_mem.ck_znQS[qmax] == null) {\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znQS[jj], cv_mem.cv_Ns);\r
+               N_VDestroyVectorArray_Serial(ck_mem.ck_znS[qmax], cv_mem.cv_Ns);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroyVectorArray_Serial(ck_mem.ck_znS[jj], cv_mem.cv_Ns);\r
+               if (ck_mem.ck_quadr) {\r
+                 N_VDestroy(ck_mem.ck_znQ[qmax]);\r
+                 for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               }\r
+               N_VDestroy(ck_mem.ck_zn[qmax]);\r
+               for (jj=0; jj<=cv_mem.cv_q; jj++) N_VDestroy(ck_mem.ck_zn[jj]);\r
+               ck_mem = null;\r
+               return(null);\r
+             }\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
+         /* Load check point data from cv_mem */\r
 \r
-                         /* Perform the actual interpolation. */\r
+         for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_zn[j], ck_mem.ck_zn[j]);\r
+         if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_zn[qmax], ck_mem.ck_zn[qmax]);\r
 \r
-                         factor1 = t - t0;\r
+         if (ck_mem.ck_quadr) {\r
+           for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znQ[j], ck_mem.ck_znQ[j]);\r
+           if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znQ[qmax], ck_mem.ck_znQ[qmax]);\r
+         }\r
 \r
-                         factor2 = factor1/delta;\r
-                         factor2 = factor2*factor2;\r
+         if (ck_mem.ck_sensi) {\r
+           for (is=0; is<cv_mem.cv_Ns; is++) {\r
+             for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znS[j][is], ck_mem.ck_znS[j][is]);\r
+             if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znS[qmax][is], ck_mem.ck_znS[qmax][is]);\r
+           }\r
+         }\r
 \r
-                         factor3 = factor2*(t-t1)/delta;\r
+         if (ck_mem.ck_quadr_sensi) {\r
+           for (is=0; is<cv_mem.cv_Ns; is++) {\r
+             for (j=0; j<=cv_mem.cv_q; j++) N_VScale_Serial(ONE, cv_mem.cv_znQS[j][is], ck_mem.ck_znQS[j][is]);\r
+             if ( cv_mem.cv_q < qmax ) N_VScale_Serial(ONE, cv_mem.cv_znQS[qmax][is], ck_mem.ck_znQS[qmax][is]);\r
+           }\r
+         }\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
+         for (j=0; j<=L_MAX; j++)        ck_mem.ck_tau[j] = cv_mem.cv_tau[j];\r
+         for (j=0; j<=NUM_TESTS; j++)    ck_mem.ck_tq[j] = cv_mem.cv_tq[j];\r
+         for (j=0; j<=cv_mem.cv_q; j++) ck_mem.ck_l[j] = cv_mem.cv_l[j];\r
+         ck_mem.ck_nst       = cv_mem.cv_nst;\r
+         ck_mem.ck_tretlast  = cv_mem.cv_tretlast;\r
+         ck_mem.ck_q         = cv_mem.cv_q;\r
+         ck_mem.ck_qprime    = cv_mem.cv_qprime;\r
+         ck_mem.ck_qwait     = cv_mem.cv_qwait;\r
+         ck_mem.ck_L         = cv_mem.cv_L;\r
+         ck_mem.ck_gammap    = cv_mem.cv_gammap;\r
+         ck_mem.ck_h         = cv_mem.cv_h;\r
+         ck_mem.ck_hprime    = cv_mem.cv_hprime;\r
+         ck_mem.ck_hscale    = cv_mem.cv_hscale;\r
+         ck_mem.ck_eta       = cv_mem.cv_eta;\r
+         ck_mem.ck_etamax    = cv_mem.cv_etamax;\r
+         ck_mem.ck_t0        = cv_mem.cv_tn;\r
+         ck_mem.ck_saved_tq5 = cv_mem.cv_saved_tq5;\r
+\r
+         return(ck_mem);\r
+       }\r
+               \r
+       /* \r
+        * CVodeGetQuad\r
+        *\r
+        * This routine extracts quadrature solution into yQout at the\r
+        * time which CVode returned the solution.\r
+        * This is just a wrapper that calls CVodeGetQuadDky with k=0.\r
+        */\r
+        \r
+       private int CVodeGetQuad(CVodeMemRec cv_mem, double tret[], NVector yQout)\r
+       {\r
+         int flag;\r
+       \r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuad", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         } \r
+       \r
+         tret[0] = cv_mem.cv_tretlast;\r
+         \r
+         flag = CVodeGetQuadDky(cv_mem,cv_mem.cv_tretlast,0,yQout);\r
+       \r
+         return(flag);\r
+       }\r
+       \r
+       /*\r
+        * CVodeGetQuadDky\r
+        *\r
+        * CVodeQuadDky computes the kth derivative of the yQ function at\r
+        * time t, where tn-hu <= t <= tn, tn denotes the current         \r
+        * internal time reached, and hu is the last internal step size   \r
+        * successfully used by the solver. The user may request \r
+        * k=0, 1, ..., qu, where qu is the current order. \r
+        * The derivative vector is returned in dky. This vector \r
+        * must be allocated by the caller. It is only legal to call this         \r
+        * function after a successful return from CVode with quadrature\r
+        * computation enabled.\r
+        */\r
+       \r
+       private int CVodeGetQuadDky(CVodeMemRec cv_mem, double t, int k, NVector dkyQ)\r
+       { \r
+         double s, c, r;\r
+         double tfuzz, tp, tn1;\r
+         int i, j;\r
+         \r
+         /* Check all inputs for legality */\r
+         \r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetQuadDky", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         } \r
+       \r
+         if(cv_mem.cv_quadr != true) {\r
+           cvProcessError(cv_mem, CV_NO_QUAD, "CVODES", "CVodeGetQuadDky", MSGCV_NO_QUAD);\r
+           return(CV_NO_QUAD);\r
+         }\r
+       \r
+         if (dkyQ == null) {\r
+           cvProcessError(cv_mem, CV_BAD_DKY, "CVODES", "CVodeGetQuadDky", MSGCV_NULL_DKY);\r
+           return(CV_BAD_DKY);\r
+         }\r
+         \r
+         if ((k < 0) || (k > cv_mem.cv_q)) {\r
+           cvProcessError(cv_mem, CV_BAD_K, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_K);\r
+           return(CV_BAD_K);\r
+         }\r
+         \r
+         /* Allow for some slack */\r
+         tfuzz = FUZZ_FACTOR * cv_mem.cv_uround *\r
+           (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_hu));\r
+         if (cv_mem.cv_hu < ZERO) tfuzz = -tfuzz;\r
+         tp = cv_mem.cv_tn - cv_mem.cv_hu - tfuzz;\r
+         tn1 = cv_mem.cv_tn + tfuzz;\r
+         if ((t-tp)*(t-tn1) > ZERO) {\r
+           cvProcessError(cv_mem, CV_BAD_T, "CVODES", "CVodeGetQuadDky", MSGCV_BAD_T);\r
+           return(CV_BAD_T);\r
+         }\r
+         \r
+         /* Sum the differentiated interpolating polynomial */\r
+         \r
+         s = (t - cv_mem.cv_tn) / cv_mem.cv_h;\r
+         for (j=cv_mem.cv_q; j >= k; j--) {\r
+           c = ONE;\r
+           for (i=j; i >= j-k+1; i--) c *= i;\r
+           if (j == cv_mem.cv_q) {\r
+             N_VScale_Serial(c, cv_mem.cv_znQ[cv_mem.cv_q], dkyQ);\r
+           } else {\r
+             N_VLinearSum_Serial(c, cv_mem.cv_znQ[j], s, dkyQ, dkyQ);\r
+           }\r
+         }\r
+         if (k == 0) return(CV_SUCCESS);\r
+         r = Math.pow(cv_mem.cv_h, -k);\r
+         N_VScale_Serial(r, dkyQ, dkyQ);\r
+         return(CV_SUCCESS);\r
+         \r
+       }\r
+       \r
+       class CVadjCheckPointRec {\r
+                 CkpntMemRec my_addr;\r
+                 CkpntMemRec next_addr;\r
+                 double t0;\r
+                 double t1;\r
+                 long nstep;\r
+                 int order;\r
+                 double step;\r
+               };\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
+        * CVodeGetAdjCheckPointsInfo\r
+        *\r
+        * This routine loads an array of nckpnts structures of type CVadjCheckPointRec.\r
+        * The user must allocate space for ckpnt.\r
+        */\r
 \r
+       private int CVodeGetAdjCheckPointsInfo(CVodeMemRec cv_mem, CVadjCheckPointRec ckpnt[])\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CkpntMemRec ck_mem;\r
+         int i;\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
-                       static 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
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeGetAdjCheckPointsInfo", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\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
+         ck_mem = ca_mem.ck_mem;\r
+\r
+         i = 0;\r
+\r
+         while (ck_mem != null) {\r
+\r
+           ckpnt[i].my_addr = ck_mem;\r
+           ckpnt[i].next_addr = ck_mem.ck_next;\r
+           ckpnt[i].t0 = ck_mem.ck_t0;\r
+           ckpnt[i].t1 = ck_mem.ck_t1;\r
+           ckpnt[i].nstep = ck_mem.ck_nst;\r
+           ckpnt[i].order = ck_mem.ck_q;\r
+           ckpnt[i].step = ck_mem.ck_h;\r
+\r
+           ck_mem = ck_mem.ck_next;\r
+           i++;\r
+\r
+         }\r
 \r
-                         if ( to_left ) {\r
-                           /* look for a new indx to the left */\r
+         return(CV_SUCCESS);\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
+       \r
+       private int CVodeCreateB(CVodeMemRec cv_mem, int lmmB, int iterB, int which[])\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec new_cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
 \r
-                           if ( indx[0] == 0 )\r
-                             ilast = 1;\r
-                           else\r
-                             ilast = indx[0];\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeCreateB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\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
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeCreateB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
 \r
-                         } else if ( to_right ) {\r
-                           /* look for a new indx to the right */\r
+         /* Allocate space for new CVodeBMem object */\r
 \r
-                           newpoint[0] = true;\r
+         new_cvB_mem = null;\r
+         new_cvB_mem = new CVodeBMemRec();\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
+         /* Create and set a new CVODES object for the backward problem */\r
 \r
-                           ilast = indx[0];\r
+         cvodeB_mem = CVodeCreate(lmmB, iterB);\r
+         if (cvodeB_mem == null) {\r
+           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeCreateB", MSGCV_MEM_FAIL);\r
+           return(CV_MEM_FAIL);\r
+         }\r
 \r
+         UserData userData = new UserData();\r
+         userData.memRec = cv_mem;\r
+         cvodeB_mem.cv_user_data = userData;\r
 \r
-                         } else {\r
-                           /* ilast is still OK */\r
+         cvodeB_mem.cv_mxhnil = -1;\r
 \r
-                           indx[0] = ilast;\r
+         //CVodeSetErrHandlerFn(cvodeB_mem, cv_mem->cv_ehfun, cv_mem->cv_eh_data);\r
+         //cvodeB_mem.cv_ehfun = cv_mem.cv_ehfun;\r
+         cvodeB_mem.cv_eh_data = cv_mem.cv_eh_data;\r
 \r
-                         }\r
+         //CVodeSetErrFile(cvodeB_mem, cv_mem->cv_errfp);\r
+         cvodeB_mem.cv_errfp = cv_mem.cv_errfp;\r
+\r
+         /* Set/initialize fields in the new CVodeBMem object, new_cvB_mem */\r
+\r
+         new_cvB_mem.cv_index   = ca_mem.ca_nbckpbs;\r
+\r
+         new_cvB_mem.cv_mem     = cvodeB_mem;\r
+\r
+         new_cvB_mem.cv_f       = -1;\r
+         //new_cvB_mem.cv_fs      = null;\r
+\r
+         new_cvB_mem.cv_fQ      = -1;\r
+         //new_cvB_mem.cv_fQs     = null;\r
+\r
+         new_cvB_mem.cv_user_data  = null;\r
+\r
+         new_cvB_mem.cv_lmem    = null;\r
+         //new_cvB_mem.cv_lfree   = null;\r
+         new_cvB_mem.cv_lfree_select = -1;\r
+         //new_cvB_mem.cv_pmem    = null;\r
+         //new_cvB_mem.cv_pfree   = null;\r
+\r
+         new_cvB_mem.cv_y       = null;\r
+\r
+         new_cvB_mem.cv_f_withSensi = false;\r
+         new_cvB_mem.cv_fQ_withSensi = false;\r
+\r
+         /* Attach the new object to the linked list cvB_mem */\r
+\r
+         new_cvB_mem.cv_next = ca_mem.cvB_mem;\r
+         ca_mem.cvB_mem = new_cvB_mem;\r
+         \r
+         /* Return the index of the newly created CVodeBMem object.\r
+          * This must be passed to CVodeInitB and to other ***B \r
+          * functions to set optional inputs for this backward problem */\r
+\r
+         which[0] = ca_mem.ca_nbckpbs;\r
+\r
+         ca_mem.ca_nbckpbs++;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       private int CVodeInitB(CVodeMemRec cv_mem, int which, \r
+               int fB,\r
+               double tB0, NVector yB0)\r
+{\r
+  CVadjMemRec ca_mem;\r
+  CVodeBMemRec cvB_mem;\r
+  CVodeMemRec cvodeB_mem;\r
+  int flag;\r
+\r
+  /* Check if cvode_mem exists */\r
+\r
+  if (cv_mem == null) {\r
+    cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeInitB", MSGCV_NO_MEM);\r
+    return(CV_MEM_NULL);\r
+  }\r
+\r
+  /* Was ASA initialized? */\r
+\r
+  if (cv_mem.cv_adjMallocDone == false) {\r
+    cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeInitB", MSGCV_NO_ADJ);\r
+    return(CV_NO_ADJ);\r
+  } \r
+  ca_mem = cv_mem.cv_adj_mem;\r
+\r
+  /* Check the value of which */\r
+\r
+  if ( which >= ca_mem.ca_nbckpbs ) {\r
+    cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeInitB", MSGCV_BAD_WHICH);\r
+    return(CV_ILL_INPUT);\r
+  }\r
+\r
+  /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+\r
+  cvB_mem = ca_mem.cvB_mem;\r
+  while (cvB_mem != null) {\r
+    if ( which == cvB_mem.cv_index ) break;\r
+    cvB_mem = cvB_mem.cv_next;\r
+  }\r
+\r
+  cvodeB_mem = cvB_mem.cv_mem;\r
+  \r
+  /* Allocate and set the CVODES object */\r
+\r
+  flag = CVodeInit(cvodeB_mem, CVArhs_select, tB0, yB0);\r
+\r
+  if (flag != CV_SUCCESS) return(flag);\r
+\r
+  /* Copy fB function in cvB_mem */\r
+\r
+  cvB_mem.cv_f_withSensi = false;\r
+  cvB_mem.cv_f = fB;\r
+\r
+  /* Allocate space and initialize the y Nvector in cvB_mem */\r
+\r
+  cvB_mem.cv_t0 = tB0;\r
+  cvB_mem.cv_y = N_VClone(yB0);\r
+  N_VScale_Serial(ONE, yB0, cvB_mem.cv_y);\r
+\r
+  return(CV_SUCCESS);\r
+}\r
+\r
+       \r
+       \r
+       \r
+\r
+       private int CVodeSStolerancesB(CVodeMemRec cv_mem, int which, double reltolB, double abstolB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeSStolerancesB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeSStolerancesB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check the value of which */\r
+\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeSStolerancesB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         cvodeB_mem = cvB_mem.cv_mem;\r
+\r
+         /* Set tolerances */\r
+\r
+         flag = CVodeSStolerances(cvodeB_mem, reltolB, abstolB);\r
+\r
+         return(flag);\r
+       }\r
+\r
+       private int CVodeSetUserDataB(CVodeMemRec cv_mem, int which, UserData user_dataB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeSetUserDataB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeSetUserDataB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeSetUserDataB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         cvB_mem.cv_user_data = user_dataB;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /*---------------------------------------------------------------\r
+         CVDlsSetLinearSolverB specifies the direct linear solver for \r
+         backward integration.\r
+         ---------------------------------------------------------------*/\r
+       private int CVDlsSetLinearSolverB(CVodeMemRec cv_mem, int which,\r
+                                 SUNLinearSolver LS, double A[][])\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         CVDlsMemRecB cvdlsB_mem;\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS",\r
+                          "CVDlsSetLinearSolverB", 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
+                          "CVDlsSetLinearSolverB", MSGD_NO_ADJ);\r
+           return(CVDLS_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS",\r
+                          "CVDlsSetLinearSolverB", MSGD_BAD_WHICH);\r
+           return(CVDLS_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         /* Get memory for CVDlsMemRecB */\r
+         cvdlsB_mem = new CVDlsMemRecB();\r
+\r
+         /* free any existing system solver attached to cvB */\r
+         //if (cvB_mem.cv_lfree_select > 0)  cvB_mem.cv_lfree(cvB_mem);\r
+\r
+         /* Attach lmemB data and lfreeB function. */\r
+         cvB_mem.cv_lmem  = cvdlsB_mem;\r
+         cvB_mem.cv_lfree_select = cvDlsFreeB_select;\r
+\r
+         /* initialize jacB and jacBS pointers */\r
+         cvdlsB_mem.jacB  = -1;\r
+         cvdlsB_mem.jacBS = -1;\r
+\r
+         /* set the linear solver for this backward problem */\r
+         cvodeB_mem = cvB_mem.cv_mem;\r
+         flag = CVDlsSetLinearSolver(cvodeB_mem, LS, A);\r
+         if (flag != CVDLS_SUCCESS) {\r
+           cvdlsB_mem = null;\r
+         }\r
+\r
+         return(flag);\r
+       }\r
+\r
+       /*-----------------------------------------------------------------\r
+         Type: CVDlsJacFnB\r
+         -----------------------------------------------------------------\r
+         A Jacobian approximation function jacB for the adjoint\r
+         (backward) problem must have the prototype given below. \r
+         -----------------------------------------------------------------*/\r
+       //typedef int (*CVDlsJacFnB)(realtype t, N_Vector y, N_Vector yB,\r
+                                  //N_Vector fyB, SUNMatrix JB,\r
+                                 // void *user_dataB, N_Vector tmp1B,\r
+                                 // N_Vector tmp2B, N_Vector tmp3B);\r
+\r
+       /*-----------------------------------------------------------------\r
+         Type: CVDlsJacFnBS\r
+         -----------------------------------------------------------------\r
+         A Jacobian approximation function jacBS for the adjoint\r
+         (backward) problem, sensitivity-dependent case,  must have the\r
+         prototype given below. \r
+         -----------------------------------------------------------------*/\r
+       //typedef int (*CVDlsJacFnBS)(realtype t, N_Vector y, N_Vector *yS,\r
+                                  // N_Vector yB, N_Vector fyB, SUNMatrix JB,\r
+                                   //void *user_dataB, N_Vector tmp1B,\r
+                                  // N_Vector tmp2B, N_Vector tmp3B);\r
+\r
+\r
+       \r
+                               \r
+       private int CVDlsSetJacFnB(CVodeMemRec cv_mem, int which, int jacB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVDlsMemRecB cvdlsB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS",\r
+                          "CVDlsSetJacFnB", 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
+                          "CVDlsSetJacFnB", MSGD_NO_ADJ);\r
+           return(CVDLS_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS",\r
+                          "CVDlsSetJacFnB", MSGD_BAD_WHICH);\r
+           return(CVDLS_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         cvodeB_mem = cvB_mem.cv_mem;\r
+\r
+         if (cvB_mem.cv_lmem == null) {\r
+           cvProcessError(cv_mem, CVDLS_LMEMB_NULL, "CVSDLS",\r
+                          "CVDlsSetJacFnB", MSGD_LMEMB_NULL);\r
+           return(CVDLS_LMEMB_NULL);\r
+         }\r
+         cvdlsB_mem = cvB_mem.cv_lmem;\r
+\r
+         cvdlsB_mem.jacB = jacB;\r
+\r
+         if (jacB > 0) {\r
+           flag = CVDlsSetJacFn(cvodeB_mem, cvDlsJacBWrapper_select);\r
+         } else {\r
+           flag = CVDlsSetJacFn(cvodeB_mem, -1);\r
+         }\r
 \r
-                         return(CV_SUCCESS);\r
+         return(flag);\r
+       }\r
+       \r
+               private int CVodeQuadInitB(CVodeMemRec cv_mem, int which,\r
+                  int fQB, NVector yQB0)\r
+       {\r
+       CVadjMemRec ca_mem;\r
+       CVodeBMemRec cvB_mem;\r
+       CVodeMemRec cvodeB_mem;\r
+       int flag;\r
+       \r
+       /* Check if cvode_mem exists */\r
+       if (cv_mem == null) {\r
+       cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeQuadInitB", MSGCV_NO_MEM);\r
+       return(CV_MEM_NULL);\r
+       }\r
+       \r
+       /* Was ASA initialized? */\r
+       if (cv_mem.cv_adjMallocDone == false) {\r
+       cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadInitB", MSGCV_NO_ADJ);\r
+       return(CV_NO_ADJ);\r
+       } \r
+       ca_mem = cv_mem.cv_adj_mem;\r
+       \r
+       /* Check which */\r
+       if ( which >= ca_mem.ca_nbckpbs ) {\r
+       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadInitB", MSGCV_BAD_WHICH);\r
+       return(CV_ILL_INPUT);\r
+       }\r
+       \r
+       /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+       cvB_mem = ca_mem.cvB_mem;\r
+       while (cvB_mem != null) {\r
+       if ( which == cvB_mem.cv_index ) break;\r
+       cvB_mem = cvB_mem.cv_next;\r
+       }\r
+       \r
+       cvodeB_mem = cvB_mem.cv_mem;\r
+       \r
+       flag = CVodeQuadInit(cvodeB_mem, CVArhsQ_select, yQB0);\r
+       if (flag != CV_SUCCESS) return(flag);\r
+       \r
+       cvB_mem.cv_fQ_withSensi = false;\r
+       cvB_mem.cv_fQ = fQB;\r
+       \r
+       return(CV_SUCCESS);\r
+       }\r
+               \r
+       /*\r
+        * CVodeSetQuad*B\r
+        *\r
+        * Wrappers for the backward phase around the corresponding \r
+        * CVODES quadrature optional input functions\r
+        */\r
+\r
+       private int CVodeSetQuadErrConB(CVodeMemRec cv_mem, int which, boolean errconQB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         int flag;\r
 \r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeSetQuadErrConB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeSetQuadErrConB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeSetQuadErrConB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         cvodeB_mem = cvB_mem.cv_mem;\r
+\r
+         cvodeB_mem.cv_errconQ = errconQB;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       private int CVodeQuadSStolerancesB(CVodeMemRec cv_mem, int which, double reltolQB, double abstolQB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem;\r
+         CVodeMemRec cvodeB_mem;\r
+         int flag;\r
+\r
+         /* Check if cvode_mem exists */\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         } \r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check which */\r
+         if ( which >= ca_mem.ca_nbckpbs ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_BAD_WHICH);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Find the CVodeBMem entry in the linked list corresponding to which */\r
+         cvB_mem = ca_mem.cvB_mem;\r
+         while (cvB_mem != null) {\r
+           if ( which == cvB_mem.cv_index ) break;\r
+           cvB_mem = cvB_mem.cv_next;\r
+         }\r
+\r
+         cvodeB_mem = cvB_mem.cv_mem;\r
+\r
+         flag = CVodeQuadSStolerances(cvodeB_mem, reltolQB, abstolQB);\r
+\r
+         return(flag);\r
+       }\r
+       \r
+       /*\r
+        * CVodeB\r
+        *\r
+        * This routine performs the backward integration towards tBout\r
+        * of all backward problems that were defined.\r
+        * When necessary, it performs a forward integration between two \r
+        * consecutive check points to update interpolation data.\r
+        *\r
+        * On a successful return, CVodeB returns CV_SUCCESS.\r
+        *\r
+        * NOTE that CVodeB DOES NOT return the solution for the backward\r
+        * problem(s). Use CVodeGetB to extract the solution at tBret\r
+        * for any given backward problem.\r
+        *\r
+        * If there are multiple backward problems and multiple check points,\r
+        * CVodeB may not succeed in getting all problems to take one step\r
+        * when called in ONE_STEP mode.\r
+        */\r
+\r
+       private int CVodeB(CVodeMemRec cv_mem, double tBout, int itaskB)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         CVodeBMemRec cvB_mem, tmp_cvB_mem;\r
+         CkpntMemRec ck_mem;\r
+         int sign, flag=0;\r
+         double tfuzz, tBn;\r
+         double tBret[] = new double[1];\r
+         boolean gotCheckpoint, isActive, reachedTBout;\r
+         \r
+         /* Check if cvode_mem exists */\r
+\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeB", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was ASA initialized? */\r
+\r
+         if (cv_mem.cv_adjMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeB", MSGCV_NO_ADJ);\r
+           return(CV_NO_ADJ);\r
+         }\r
+         ca_mem = cv_mem.cv_adj_mem;\r
+\r
+         /* Check if any backward problem has been defined */\r
+\r
+         if ( ca_mem.ca_nbckpbs == 0 ) {\r
+           cvProcessError(cv_mem, CV_NO_BCK, "CVODEA", "CVodeB", MSGCV_NO_BCK);\r
+           return(CV_NO_BCK);\r
+         }\r
+         cvB_mem = ca_mem.cvB_mem;\r
+\r
+         /* Check whether CVodeF has been called */\r
+\r
+         if ( ca_mem.ca_firstCVodeFcall ) {\r
+           cvProcessError(cv_mem, CV_NO_FWD, "CVODEA", "CVodeB", MSGCV_NO_FWD);\r
+           return(CV_NO_FWD);\r
+         }\r
+         sign = (ca_mem.ca_tfinal - ca_mem.ca_tinitial > ZERO) ? 1 : -1;\r
+\r
+         /* If this is the first call, loop over all backward problems and\r
+          *   - check that tB0 is valid\r
+          *   - check that tBout is ahead of tB0 in the backward direction\r
+          *   - check whether we need to interpolate forward sensitivities\r
+          */\r
+\r
+         if ( ca_mem.ca_firstCVodeBcall ) {\r
+\r
+           tmp_cvB_mem = cvB_mem;\r
+\r
+           while(tmp_cvB_mem != null) {\r
+\r
+             tBn = tmp_cvB_mem.cv_mem.cv_tn;\r
+\r
+             if ( (sign*(tBn-ca_mem.ca_tinitial) < ZERO) || (sign*(ca_mem.ca_tfinal-tBn) < ZERO) ) {\r
+               cvProcessError(cv_mem, CV_BAD_TB0, "CVODEA", "CVodeB", MSGCV_BAD_TB0,\r
+                              tmp_cvB_mem.cv_index);\r
+               return(CV_BAD_TB0);\r
+             }\r
+\r
+             if (sign*(tBn-tBout) <= ZERO) {\r
+               cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_TBOUT,\r
+                              tmp_cvB_mem.cv_index);\r
+               return(CV_ILL_INPUT);\r
+             }\r
+\r
+             if ( tmp_cvB_mem.cv_f_withSensi || tmp_cvB_mem.cv_fQ_withSensi )\r
+                 ca_mem.ca_IMinterpSensi = true;\r
+\r
+             tmp_cvB_mem = tmp_cvB_mem.cv_next;\r
+\r
+           }\r
+\r
+           if ( ca_mem.ca_IMinterpSensi && !ca_mem.ca_IMstoreSensi) {\r
+             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_SENSI);\r
+             return(CV_ILL_INPUT);\r
+           }\r
+\r
+           ca_mem.ca_firstCVodeBcall = false;\r
+         }\r
 \r
-                       }\r
+         /* Check if itaskB is legal */\r
 \r
+         if ( (itaskB != CV_NORMAL) && (itaskB != CV_ONE_STEP) ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_ITASKB);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Check if tBout is legal */\r
+\r
+         if ( (sign*(tBout-ca_mem.ca_tinitial) < ZERO) || (sign*(ca_mem.ca_tfinal-tBout) < ZERO) ) {\r
+           tfuzz = HUNDRED*cv_mem.cv_uround*(Math.abs(ca_mem.ca_tinitial) + Math.abs(ca_mem.ca_tfinal));\r
+           if ( (sign*(tBout-ca_mem.ca_tinitial) < ZERO) && (Math.abs(tBout-ca_mem.ca_tinitial) < tfuzz) ) {\r
+             tBout = ca_mem.ca_tinitial;\r
+           } else {\r
+             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_TBOUT);\r
+             return(CV_ILL_INPUT);\r
+           }\r
+         }\r
+\r
+         /* Loop through the check points and stop as soon as a backward\r
+          * problem has its tn value behind the current check point's t0_\r
+          * value (in the backward direction) */\r
+\r
+         ck_mem = ca_mem.ck_mem;\r
+\r
+         gotCheckpoint = false;\r
+\r
+         for(;;) {\r
+\r
+           tmp_cvB_mem = cvB_mem;\r
+           while(tmp_cvB_mem != null) {\r
+             tBn = tmp_cvB_mem.cv_mem.cv_tn;\r
+\r
+             if ( sign*(tBn-ck_mem.ck_t0) > ZERO ) {\r
+               gotCheckpoint = true;\r
+               break;\r
+             }\r
+\r
+             if ( (itaskB==CV_NORMAL) && (tBn == ck_mem.ck_t0) && (sign*(tBout-ck_mem.ck_t0) >= ZERO) ) {\r
+               gotCheckpoint = true;\r
+               break;\r
+             }\r
+\r
+             tmp_cvB_mem = tmp_cvB_mem.cv_next;\r
+           }\r
+\r
+           if (gotCheckpoint) break;\r
+\r
+           if (ck_mem.ck_next == null) break;\r
+\r
+           ck_mem = ck_mem.ck_next;\r
+         }\r
+\r
+         /* Starting with the current check point from above, loop over check points\r
+            while propagating backward problems */\r
+\r
+         for(;;) {\r
+\r
+           /* Store interpolation data if not available.\r
+              This is the 2nd forward integration pass */\r
+\r
+           if (ck_mem != ca_mem.ca_ckpntData) {\r
+             flag = CVAdataStore(cv_mem, ck_mem);\r
+             if (flag != CV_SUCCESS) break;\r
+           }\r
+\r
+           /* Loop through all backward problems and, if needed,\r
+            * propagate their solution towards tBout */\r
+\r
+           tmp_cvB_mem = cvB_mem;\r
+           while (tmp_cvB_mem != null) {\r
+\r
+             /* Decide if current backward problem is "active" in this check point */\r
+\r
+             isActive = true;\r
+\r
+             tBn = tmp_cvB_mem.cv_mem.cv_tn;\r
+\r
+             if ( (tBn == ck_mem.ck_t0) && (sign*(tBout-ck_mem.ck_t0) < ZERO ) ) isActive = false;\r
+             if ( (tBn == ck_mem.ck_t0) && (itaskB==CV_ONE_STEP) ) isActive = false;\r
+\r
+             if ( sign * (tBn - ck_mem.ck_t0) < ZERO ) isActive = false;\r
+\r
+             if ( isActive ) {\r
+\r
+               /* Store the address of current backward problem memory \r
+                * in ca_mem to be used in the wrapper functions */\r
+               ca_mem.ca_bckpbCrt = tmp_cvB_mem;\r
+\r
+               /* Integrate current backward problem */\r
+               CVodeSetStopTime(tmp_cvB_mem.cv_mem, ck_mem.ck_t0);\r
+               flag = CVode(tmp_cvB_mem.cv_mem, tBout, tmp_cvB_mem.cv_y, tBret, itaskB);\r
+\r
+               /* Set the time at which we will report solution and/or quadratures */\r
+               tmp_cvB_mem.cv_tout = tBret[0];\r
+\r
+               /* If an error occurred, exit while loop */\r
+               if (flag < 0) break;\r
+\r
+             } else {\r
+               flag = CV_SUCCESS;\r
+               tmp_cvB_mem.cv_tout = tBn;\r
+             }\r
+\r
+             /* Move to next backward problem */\r
+\r
+             tmp_cvB_mem = tmp_cvB_mem.cv_next;\r
+           }\r
+\r
+           /* If an error occurred, return now */\r
+\r
+           if (flag <0) {\r
+             cvProcessError(cv_mem, flag, "CVODEA", "CVodeB", MSGCV_BACK_ERROR,\r
+                            tmp_cvB_mem.cv_index);\r
+             return(flag);\r
+           }\r
+\r
+           /* If in CV_ONE_STEP mode, return now (flag = CV_SUCCESS) */\r
+\r
+           if (itaskB == CV_ONE_STEP) break;\r
+\r
+           /* If all backward problems have succesfully reached tBout, return now */\r
+\r
+           reachedTBout = true;\r
+\r
+           tmp_cvB_mem = cvB_mem;\r
+           while(tmp_cvB_mem != null) {\r
+             if ( sign*(tmp_cvB_mem.cv_tout - tBout) > ZERO ) {\r
+               reachedTBout = false;\r
+               break;\r
+             }\r
+             tmp_cvB_mem = tmp_cvB_mem.cv_next;\r
+           }\r
+\r
+           if ( reachedTBout ) break;\r
+\r
+           /* Move check point in linked list to next one */\r
+\r
+           ck_mem = ck_mem.ck_next;\r
+\r
+         } \r
+\r
+         return(flag);\r
+       }\r
+\r
+       /*\r
+        * CVAdataStore\r
+        *\r
+        * This routine integrates the forward model starting at the check\r
+        * point ck_mem and stores y and yprime at all intermediate steps.\r
+        *\r
+        * Return values:\r
+        * CV_SUCCESS\r
+        * CV_REIFWD_FAIL\r
+        * CV_FWD_FAIL\r
+        */\r
+\r
+       private int CVAdataStore(CVodeMemRec cv_mem, CkpntMemRec ck_mem)\r
+       {\r
+      CVadjMemRec ca_mem;\r
+         DtpntMemRec dt_mem[];\r
+         double t[] = new double[1];\r
+         long i;\r
+         int flag, sign;\r
+\r
+         ca_mem = cv_mem.cv_adj_mem;\r
+         dt_mem = ca_mem.dt_mem;\r
+\r
+         /* Initialize cv_mem with data from ck_mem */\r
+         //flag = CVAckpntGet(cv_mem, ck_mem);\r
+         //if (flag != CV_SUCCESS)\r
+           //return(CV_REIFWD_FAIL);\r
+\r
+         /* Set first structure in dt_mem[0] */\r
+         dt_mem[0].t = ck_mem.ck_t0;\r
+         //ca_mem.ca_IMstore(cv_mem, dt_mem[0]);\r
+\r
+         /* Decide whether TSTOP must be activated */\r
+         if (ca_mem.ca_tstopCVodeFcall) {\r
+           CVodeSetStopTime(cv_mem, ca_mem.ca_tstopCVodeF);\r
+         }\r
+\r
+         sign = (ca_mem.ca_tfinal - ca_mem.ca_tinitial > ZERO) ? 1 : -1;\r
+\r
+\r
+         /* Run CVode to set following structures in dt_mem[i] */\r
+         i = 1;\r
+         do {\r
+\r
+           flag = CVode(cv_mem, ck_mem.ck_t1, ca_mem.ca_ytmp, t, CV_ONE_STEP);\r
+           if (flag < 0) return(CV_FWD_FAIL);\r
+\r
+           dt_mem[(int)i].t = t[0];\r
+           //ca_mem.ca_IMstore(cv_mem, dt_mem[(int)i]);\r
+           i++;\r
+\r
+         } while ( sign*(ck_mem.ck_t1 - t[0]) > ZERO );\r
+\r
+\r
+         ca_mem.ca_IMnewData = true;     /* New data is now available    */\r
+         ca_mem.ca_ckpntData = ck_mem;   /* starting at this check point */\r
+         ca_mem.ca_np = i;               /* and we have this many points */\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /* \r
+        * CVodeSetStopTime\r
+        *\r
+        * Specifies the time beyond which the integration is not to proceed.\r
+        */\r
+\r
+       private int CVodeSetStopTime(CVodeMemRec cv_mem, double tstop)\r
+       {\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSetStopTime", MSGCV_NO_MEM);\r
+           return (CV_MEM_NULL);\r
+         }\r
+\r
+         /* If CVode was called at least once, test if tstop is legal\r
+          * (i.e. if it was not already passed).\r
+          * If CVodeSetStopTime is called before the first call to CVode,\r
+          * tstop will be checked in CVode. */\r
+         if (cv_mem.cv_nst > 0) {\r
+\r
+           if ( (tstop - cv_mem.cv_tn) * cv_mem.cv_h < ZERO ) {\r
+             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSetStopTime", MSGCV_BAD_TSTOP, tstop, cv_mem.cv_tn);\r
+             return(CV_ILL_INPUT);\r
+           }\r
+\r
+         }\r
+\r
+         cv_mem.cv_tstop = tstop;\r
+         cv_mem.cv_tstopset = true;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
 \r
 \r
 \r