File becoming too large. Must put runCVSRoberts_ASAi_dns() in a second file if work...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 2 Apr 2018 21:27:26 +0000 (21:27 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 2 Apr 2018 21:27:26 +0000 (21:27 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15434 ba61647d-9d00-f842-95cd-605cb4296b96

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

index c8fb0c6..da53192 100644 (file)
@@ -123,6 +123,21 @@ public abstract class CVODES {
        final int CV_STAGGERED1 = 3;\r
     int sensi_meth = CV_SIMULTANEOUS;\r
     boolean err_con = true;\r
+    \r
+    /* \r
+     * ----------------------------------------\r
+     * CVODEA return flags\r
+     * ----------------------------------------\r
+     */\r
+\r
+    final int CV_NO_ADJ = -101;\r
+    final int CV_NO_FWD = -102;\r
+    final int CV_NO_BCK = -103;\r
+    final int CV_BAD_TB0 = -104;\r
+    final int CV_REIFWD_FAIL = -105;\r
+    final int CV_FWD_FAIL = -106;\r
+    final int CV_GETY_BADT = -107;\r
+\r
 \r
        \r
        /*\r
@@ -155,6 +170,10 @@ public abstract class CVODES {
        /* interp */\r
        final int CV_HERMITE = 1;\r
        final int CV_POLYNOMIAL = 2;\r
+       final int CVAhermiteMalloc_select = 1;\r
+       final int CVApolynomialMalloc_select = 2;\r
+       final int CVAhermiteStorePnt_select = 1;\r
+       final int CVApolynomialStorePnt_select = 2;\r
 \r
        /* \r
         * ----------------------------------------\r
@@ -689,9 +708,10 @@ public abstract class CVODES {
            else if (problem == cvsRoberts_FSA_dns) {\r
                runcvsRoberts_FSA_dns();\r
            }\r
-           else if (problem == cvsRoberts_ASAi_dns) {\r
-               runcvsRoberts_ASAi_dns();\r
-           }\r
+           //else if (problem == cvsRoberts_ASAi_dns) {\r
+               //runcvsRoberts_ASAi_dns();\r
+               // File becoming too large.  Must put runcvsRoberts_ASAi_dns() in second file to continue.\r
+           //}\r
            else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                runcvsDirectDemo_Problem_1();\r
            }\r
@@ -1851,6 +1871,8 @@ public abstract class CVODES {
                double A[][];\r
                SUNLinearSolver LS;\r
                int steps;\r
+               double time[] = new double[1];\r
+               int ncheck[] = new int[1];\r
                \r
                // Print problem description \r
                System.out.printf("\nAdjoint Sensitivity Example for Chemical Kinetics\n");\r
@@ -1973,14 +1995,26 @@ public abstract class CVODES {
                return;\r
         }\r
 \r
-        /* Call CVodeAdjInit to update CVODES memory block by allocting the internal \r
+        /* Call CVodeAdjInit to update CVODES memory block by allocating the internal \r
         memory needed for backward integration.*/\r
-        steps = STEPS; /* no. of integration steps between two consecutive ckeckpoints*/\r
+        steps = STEPS; /* no. of integration steps between two consecutive checkpoints*/\r
         flag = CVodeAdjInit(cvode_mem, steps, CV_HERMITE);\r
         /*\r
         flag = CVodeAdjInit(cvode_mem, steps, CV_POLYNOMIAL);\r
         */\r
-\r
+        if (flag != CV_SUCCESS) {\r
+               return;\r
+        }\r
+        \r
+        // Perform forward run\r
+        System.out.printf("Forward integration ... ");\r
+       \r
+        // Call CVodeF to integrate the forward problem over an interval in time and\r
+        // save checkpointing data.\r
+        flag = CVodeF(cvode_mem, TOUT, y, time, CV_NORMAL, ncheck);\r
+        if (flag < 0) {\r
+               return;\r
+        }\r
 \r
        }\r
 \r
@@ -12205,10 +12239,6 @@ else                return(snrm);
 \r
          ca_mem = null;\r
          ca_mem = new CVadjMemRec();\r
-         if (ca_mem == null) {\r
-           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);\r
-           return(CV_MEM_FAIL);\r
-         }\r
 \r
          /* Attach ca_mem to CVodeMem structure */\r
 \r
@@ -12267,19 +12297,19 @@ else                return(snrm);
 \r
          case CV_HERMITE:\r
            \r
-           //ca_mem.ca_IMmalloc = CVAhermiteMalloc;\r
+           ca_mem.ca_IMmalloc = CVAhermiteMalloc_select;\r
            //ca_mem.ca_IMfree   = CVAhermiteFree;\r
            //ca_mem.ca_IMget    = CVAhermiteGetY;\r
-           //ca_mem.ca_IMstore  = CVAhermiteStorePnt;\r
+           ca_mem.ca_IMstore  = CVAhermiteStorePnt_select;\r
 \r
            break;\r
            \r
          case CV_POLYNOMIAL:\r
          \r
-           //ca_mem.ca_IMmalloc = CVApolynomialMalloc;\r
+           ca_mem.ca_IMmalloc = CVApolynomialMalloc_select;\r
            //ca_mem.ca_IMfree   = CVApolynomialFree;\r
            //ca_mem.ca_IMget    = CVApolynomialGetY;\r
-           //ca_mem.ca_IMstore  = CVApolynomialStorePnt;\r
+           ca_mem.ca_IMstore  = CVApolynomialStorePnt_select;\r
 \r
            break;\r
 \r
@@ -12387,9 +12417,9 @@ else                return(snrm);
                  int ca_IMtype;\r
 \r
                  /* Functions set by the interpolation module */\r
-                 //cvaIMMallocFn   ca_IMmalloc; \r
+                 int   ca_IMmalloc; \r
                  //cvaIMFreeFn     ca_IMfree;\r
-                 //cvaIMStorePntFn ca_IMstore; /* store a new interpolation point */\r
+                 int ca_IMstore; /* store a new interpolation point */\r
                  //cvaIMGetYFn     ca_IMget;   /* interpolate forward solution    */\r
 \r
                  /* Flags controlling the interpolation module */\r
@@ -12555,9 +12585,848 @@ else                return(snrm);
                  \r
                class DtpntMemRec {\r
                  double t;    /* time */\r
-                 //void *content; /* IMtype-dependent content */\r
+                 int content;\r
+                 HermiteDataMemRec hermiteContent;\r
+                 PolynomialDataMemRec polynomialContent;\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
+                 /* 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
+                 /* 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
+\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
+\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
+\r
+                 /* ck_mem->ck_zn[qmax] was not allocated */\r
+                 ck_mem.ck_zqm = 0;\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
+                 if (ck_mem.ck_quadr) {\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
+                   N_VScale_Serial(ONE, cv_mem.cv_znQ[0], ck_mem.ck_znQ[0]);\r
+\r
+                 }\r
+\r
+                 /* Do 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
+                   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
+                   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
+                 /* 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
+\r
+                 /* Next in list */\r
+                 ck_mem.ck_next  = null;\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
+\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
+                 /* 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
+\r
+                 }\r
+\r
+                 return(allocOK);\r
+               }\r
+\r
+               /* Data for cubic Hermite interpolation */\r
+               class HermiteDataMemRec {\r
+                 NVector y;\r
+                 NVector yd;\r
+                 NVector yS[];\r
+                 NVector ySd[];\r
+               };\r
+\r
+               /* Data for polynomial interpolation */\r
+               class PolynomialDataMemRec {\r
+                 NVector y;\r
+                 NVector yS[];\r
+                 int order;\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
+\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
+                 /* 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 PolynomialDataMemRec();\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
+                   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
+                       content = null;\r
+                       ii = i;\r
+                       allocOK = false;\r
+                       break;\r
+                     }\r
+\r
+                   }\r
+\r
+                   dt_mem[i].content = CV_POLYNOMIAL;\r
+                   dt_mem[i].polynomialContent = 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].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
+                 }\r
+\r
+                 return(allocOK);\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
+\r
+                 /* Load derivative */\r
+\r
+                 if (cv_mem.cv_nst == 0) {\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
+\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
+                 } else {\r
+\r
+                   N_VScale_Serial(ONE/cv_mem.cv_h, cv_mem.cv_zn[1], content.yd);\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
+\r
+                 }\r
+\r
+                 return(0);\r
+               }\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
+               private int CVApolynomialStorePnt(CVodeMemRec cv_mem, DtpntMemRec d)\r
+               {\r
+                 CVadjMemRec ca_mem;\r
+                 PolynomialDataMemRec content;\r
+                 int is;\r
+\r
+                 ca_mem = cv_mem.cv_adj_mem;\r
+\r
+                 content = d.polynomialContent;\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
+                 content.order = cv_mem.cv_qu;\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
+                 if (ck_mem == null) return(null);\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 (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
+                 /* 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
+\r
+                 }\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
+\r
+                 }\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
+\r
+                 }\r
+\r
+                 /* Load check point data from cv_mem */\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
+\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
+                 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
+                 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
+                 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
 \r
 }
\ No newline at end of file