Created for running runcvsRoberts_ASAi_dns().
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 18:25:16 +0000 (18:25 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 18:25:16 +0000 (18:25 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15436 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/CVODES_ASA.java [new file with mode: 0644]

diff --git a/mipav/src/gov/nih/mipav/model/algorithms/CVODES_ASA.java b/mipav/src/gov/nih/mipav/model/algorithms/CVODES_ASA.java
new file mode 100644 (file)
index 0000000..0ef94ce
--- /dev/null
@@ -0,0 +1,1437 @@
+package gov.nih.mipav.model.algorithms;\r
+\r
+\r
+import gov.nih.mipav.view.MipavUtil;\r
+\r
+public abstract class CVODES_ASA extends CVODES {\r
+       public CVODES_ASA() {\r
+               super();\r
+               runcvsRoberts_ASAi_dns();\r
+       }\r
+       \r
+       // In CVODES.java have int problem = cvsRoberts_ASAi_dns;\r
+       // Use the following code to call CVODES_ASA from another module:\r
+               /*boolean testme = true;\r
+           class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns()\r
+                 public CVODEStest() {\r
+                         super();\r
+                 }\r
+                 public int f(double t, NVector yv, NVector ydotv, UserData user_data) {\r
+                         return 0;\r
+                 }\r
+                 \r
+                 public int g(double t, NVector yv, double gout[], UserData user_data) {\r
+                         return 0;\r
+                 }\r
+                 \r
+                 public int fQ(double t, NVector x, NVector y, UserData user_data) {\r
+                         return 0;\r
+                 }\r
+                 \r
+                 public int fS1(int Ns, double t, NVector yv, NVector ydot, int is,\r
+                                        NVector yS, NVector ySdot, UserData user_data, NVector tmp1, NVector tmp2) {\r
+                                 return 0;  \r
+                 }\r
+                 \r
+                 public int Jac(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, \r
+                                 NVector tmp2, NVector tmp3) {\r
+                         return 0;\r
+                 }\r
+                 \r
+                 public int ewt(NVector y, NVector w, UserData user_data) {\r
+                     return 0;\r
+                 }\r
+                 \r
+                 public int fB(double t, NVector y, NVector yB, NVector yBdot, UserData user_dataB) {\r
+                     return 0;\r
+                 }\r
+                 \r
+                 public int JacB(double t, NVector y, NVector yB, NVector fyB, double JB[][], UserData user_dataB,\r
+                               NVector tmp1B, NVector tmp2B, NVector tmp3B) {\r
+                                 return 0;\r
+                         }\r
+                 \r
+                 public int fQB(double t, NVector y, NVector yB, NVector qBdot, UserData user_dataB) {\r
+                     return 0;\r
+                 }\r
+             }\r
+           if (testme) {\r
+               new CVODEStest();\r
+               return;\r
+           } */\r
+       \r
+       private void runcvsRoberts_ASAi_dns() {\r
+               /* -----------------------------------------------------------------\r
+                * Programmer(s): Radu Serban @ LLNL\r
+                * -----------------------------------------------------------------\r
+                * LLNS Copyright Start\r
+                * Copyright (c) 2014, Lawrence Livermore National Security\r
+                * This work was performed under the auspices of the U.S. Department\r
+                * of Energy by Lawrence Livermore National Laboratory in part under\r
+                * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.\r
+                * Produced at the Lawrence Livermore National Laboratory.\r
+                * All rights reserved.\r
+                * For details, see the LICENSE file.\r
+                * LLNS Copyright End\r
+                * -----------------------------------------------------------------\r
+                * Adjoint sensitivity example problem.\r
+                * The following is a simple example problem, with the coding\r
+                * needed for its solution by CVODES. The problem is from chemical\r
+                * kinetics, and consists of the following three rate equations.\r
+                *    dy1/dt = -p1*y1 + p2*y2*y3\r
+                *    dy2/dt =  p1*y1 - p2*y2*y3 - p3*(y2)^2\r
+                *    dy3/dt =  p3*(y2)^2\r
+                * on the interval from t = 0.0 to t = 4.e10, with initial\r
+                * conditions: y1 = 1.0, y2 = y3 = 0. The reaction rates are:\r
+                * p1=0.04, p2=1e4, and p3=3e7. The problem is stiff.\r
+                * This program solves the problem with the BDF method, Newton\r
+                * iteration with the CVODE dense linear solver, and a user-supplied\r
+                * Jacobian routine.\r
+                * It uses a scalar relative tolerance and a vector absolute\r
+                * tolerance.\r
+                * Output is printed in decades from t = .4 to t = 4.e10.\r
+                * Run statistics (optional outputs) are printed at the end.\r
+                * \r
+                * Optionally, CVODES can compute sensitivities with respect to\r
+                * the problem parameters p1, p2, and p3 of the following quantity:\r
+                *   G = int_t0^t1 g(t,p,y) dt\r
+                * where\r
+                *   g(t,p,y) = y3\r
+                *        \r
+                * The gradient dG/dp is obtained as:\r
+                *   dG/dp = int_t0^t1 (g_p - lambda^T f_p ) dt - lambda^T(t0)*y0_p\r
+                *         = - xi^T(t0) - lambda^T(t0)*y0_p\r
+                * where lambda and xi are solutions of:\r
+                *   d(lambda)/dt = - (f_y)^T * lambda - (g_y)^T\r
+                *   lambda(t1) = 0\r
+                * and\r
+                *   d(xi)/dt = - (f_p)^T * lambda + (g_p)^T\r
+                *   xi(t1) = 0\r
+                * \r
+                * During the backward integration, CVODES also evaluates G as\r
+                *   G = - phi(t0)\r
+                * where\r
+                *   d(phi)/dt = g(t,y,p)\r
+                *   phi(t1) = 0\r
+                * -----------------------------------------------------------------*/\r
+       \r
+               /** Problem Constants */\r
+               final int NEQ = 3; // Number of equations\r
+               //final double RTOL = 1.0E-6; // scalar relative tolerance\r
+               final double RTOL = 1.0E-12;\r
+               //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+               final double ATOL1 = 1.0E-12;\r
+               //final double ATOL2 = 1.0E-14;\r
+               final double ATOL2 = 1.0E-15;\r
+               //final double ATOL3 = 1.0E-6;\r
+               final double ATOL3 = 1.0E-12;\r
+               final double ATOLl = 1.0E-8; // absolute tolerance for adjoint variables\r
+               final double ATOLq = 1.0E-6; // absolute tolerannce for quuadratures\r
+               final double T0 = 0.0; // initial time\r
+               final double TOUT = 4.0E7; // final time\r
+               final double TB1 = 4.0E7; // starting point for adjoint problem\r
+               final double TB2 = 50.0; // starting point for adjoint problem\r
+               final double TBout1 = 40.0; // intermediate t for adjoint problem\r
+               final int STEPS = 150; // number of steps between check points\r
+               final int NP = 3; // number of problem parameters\r
+               double reltolQ;\r
+               double abstolQ;\r
+               CVodeMemRec cvode_mem;\r
+               int flag;\r
+        int f = cvsRoberts_ASAi_dns;\r
+        int Jac = cvsRoberts_ASAi_dns;\r
+               int ewt_select = cvEwtUser_select1;\r
+               int fQ= cvsRoberts_ASAi_dns;\r
+               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
+               System.out.printf("-------------------------------------------------\n\n");\r
+               System.out.printf("ODE: dy1/dt = -p1*y1 + p2*y2*y3\n");\r
+               System.out.printf("     dy2/dt =  p1*y1 - p2*y2*y3 - p3*(y2)^2\n");\r
+               System.out.printf("     dy3/dt =  p3*(y2)^2\n\n");\r
+               System.out.printf("Find dG/dp for\n");\r
+               System.out.printf("     G = int_t0^tB0 g(t,p,y) dt\n");\r
+               System.out.printf("     g(t,p,y) = y3\n\n\n");\r
+\r
+               // User data structure\r
+               UserData data = new UserData();\r
+               data.array = new double[3];\r
+               data.array[0] = 0.04;\r
+               data.array[1] = 1.0E4;\r
+               data.array[2] = 3.0E7;\r
+               \r
+               // Initialize y\r
+               NVector y = new NVector();\r
+               double yr[] = new double[]{1.0,0.0,0.0};\r
+               N_VNew_Serial(y, NEQ);\r
+               y.setData(yr);\r
+               \r
+               // Initialize q\r
+               NVector q = new NVector();\r
+               double qr[] = new double[]{0.0};\r
+               N_VNew_Serial(q,1);\r
+               q.setData(qr);\r
+               \r
+               // Set the scalar relative and absolute tolerances reltolQ and abstolQ;\r
+               reltolQ = RTOL;\r
+               abstolQ = ATOLq;\r
+               \r
+               /* Create and allocate CVODES memory for forward run */\r
+               System.out.printf("Create and allocate CVODES memory for forward runs\n");\r
+\r
+               /* Call CVodeCreate to create the solver memory and specify the \r
+                    Backward Differentiation Formula and the use of a Newton iteration */\r
+               cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);\r
+               if (cvode_mem == null) {\r
+                   return;     \r
+               }\r
+               // Allow unlimited steps in reaching tout\r
+               flag = CVodeSetMaxNumSteps(cvode_mem, -1);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               // Allow many error test failures\r
+               flag = CVodeSetMaxErrTestFails(cvode_mem, Integer.MAX_VALUE);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeInit to initialize the integrator memory and specify the\r
+               // user's right hand side function in y' = f(t,y), the initial time T0, and\r
+           // the initial dependent variable vector y\r
+               flag = CVodeInit(cvode_mem, f, T0, y);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Use private function to compute error weights\r
+               // CVodeWFtolerances specifies a user-provides function (of type CVEwtFn)\r
+               // which will be called to set the error weight vector.\r
+               flag = CVodeWFtolerances(cvode_mem, ewt_select);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Attach user data\r
+               cvode_mem.cv_user_data = data;\r
+               \r
+               \r
+               // Create dense SUNMATRIX for use in linear solver\r
+               // indexed by columns stacked on top of each other\r
+               try {\r
+                   A = new double[NEQ][NEQ];\r
+               }\r
+               catch (OutOfMemoryError e) {\r
+                   MipavUtil.displayError("Out of memory error trying to create A");\r
+                   return;\r
+               }\r
+               \r
+               // Create dense SUNLinearSolver object for use by CVode\r
+               LS = SUNDenseLinearSolver(y, A);\r
+               if (LS == null) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
+               flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the user-supplied Jacobian routine Jac\r
+               flag = CVDlsSetJacFn(cvode_mem, Jac);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeQuadInit to allocate internal memory and initialize\r
+               // quadrature integration\r
+               flag = CVodeQuadInit(cvode_mem, fQ, q);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               /* Call CVodeSetQuadErrCon to specify whether or not the quadrature variables\r
+            are to be used in the step size control mechanism within CVODES. Call\r
+            CVodeQuadSStolerances or CVodeQuadSVtolerances to specify the integration\r
+            tolerances for the quadrature variables. */\r
+               cvode_mem.cv_errconQ = true;\r
+               \r
+               /* Call CVodeQuadSStolerances to specify scalar relative and absolute\r
+            tolerances. */\r
+           flag = CVodeQuadSStolerances(cvode_mem, reltolQ, abstolQ);\r
+        if (flag != CV_SUCCESS) {\r
+               return;\r
+        }\r
+\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 checkpoints*/\r
+        flag = CVodeAdjInit(cvode_mem, steps, CV_HERMITE);\r
+        /*\r
+        flag = CVodeAdjInit(cvode_mem, steps, CV_POLYNOMIAL);\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
+        * CVodeQuadInit\r
+        *\r
+        * CVodeQuadInit allocates and initializes quadrature related \r
+        * memory for a problem. All problem specification inputs are \r
+        * checked for errors. If any error occurs during initialization, \r
+        * it is reported to the file whose file pointer is errfp. \r
+        * The return value is CV_SUCCESS = 0 if no errors occurred, or\r
+        * a negative value otherwise.\r
+        */\r
+\r
+       private int CVodeQuadInit(CVodeMemRec cv_mem, int fQ, NVector yQ0)\r
+       {\r
+         boolean allocOK;\r
+         long lrw1Q[] = new long[1];\r
+         long liw1Q[] = new long[1];\r
+\r
+         /* Check cvode_mem */\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeQuadInit",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Set space requirements for one N_Vector */\r
+         N_VSpace_Serial(yQ0, lrw1Q, liw1Q);\r
+         cv_mem.cv_lrw1Q = lrw1Q[0];\r
+         cv_mem.cv_liw1Q = liw1Q[0];\r
+\r
+         /* Allocate the vectors (using yQ0 as a template) */\r
+         allocOK = cvQuadAllocVectors(cv_mem, yQ0);\r
+         if (!allocOK) {\r
+           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES",\r
+                          "CVodeQuadInit", MSGCV_MEM_FAIL);\r
+           return(CV_MEM_FAIL);\r
+         }\r
+\r
+         /* Initialize znQ[0] in the history array */\r
+         N_VScale_Serial(ONE, yQ0, cv_mem.cv_znQ[0]);\r
+\r
+         /* Copy the input parameters into CVODES state */\r
+         cv_mem.cv_fQ = fQ;\r
+\r
+         /* Initialize counters */\r
+         cv_mem.cv_nfQe  = 0;\r
+         cv_mem.cv_netfQ[0] = 0;\r
+\r
+         /* Quadrature integration turned ON */\r
+         cv_mem.cv_quadr = true;\r
+         cv_mem.cv_QuadMallocDone = true;\r
+\r
+         /* Quadrature initialization was successfull */\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /*\r
+        * CVodeQuadAllocVectors\r
+        *\r
+        * NOTE: Space for ewtQ is allocated even when errconQ=SUNFALSE, \r
+        * although in this case, ewtQ is never used. The reason for this\r
+        * decision is to allow the user to re-initialize the quadrature\r
+        * computation with errconQ=SUNTRUE, after an initialization with\r
+        * errconQ=SUNFALSE, without new memory allocation within \r
+        * CVodeQuadReInit.\r
+        */\r
+\r
+       private boolean cvQuadAllocVectors(CVodeMemRec cv_mem, NVector tmpl) \r
+       {\r
+         int i, j;\r
+\r
+         /* Allocate ewtQ */\r
+         cv_mem.cv_ewtQ = N_VClone(tmpl);\r
+         if (cv_mem.cv_ewtQ == null) {\r
+           return(false);\r
+         }\r
+         \r
+         /* Allocate acorQ */\r
+         cv_mem.cv_acorQ = N_VClone(tmpl);\r
+         if (cv_mem.cv_acorQ == null) {\r
+           N_VDestroy(cv_mem.cv_ewtQ);\r
+           return(false);\r
+         }\r
+\r
+         /* Allocate yQ */\r
+         cv_mem.cv_yQ = N_VClone(tmpl);\r
+         if (cv_mem.cv_yQ == null) {\r
+           N_VDestroy(cv_mem.cv_ewtQ);\r
+           N_VDestroy(cv_mem.cv_acorQ);\r
+           return(false);\r
+         }\r
+\r
+         /* Allocate tempvQ */\r
+         cv_mem.cv_tempvQ = N_VClone(tmpl);\r
+         if (cv_mem.cv_tempvQ == null) {\r
+           N_VDestroy(cv_mem.cv_ewtQ);\r
+           N_VDestroy(cv_mem.cv_acorQ);\r
+           N_VDestroy(cv_mem.cv_yQ);\r
+           return(false);\r
+         }\r
+\r
+         /* Allocate zQn[0] ... zQn[maxord] */\r
+\r
+         for (j=0; j <= cv_mem.cv_qmax; j++) {\r
+           cv_mem.cv_znQ[j] = N_VClone(tmpl);\r
+           if (cv_mem.cv_znQ[j] == null) {\r
+             N_VDestroy(cv_mem.cv_ewtQ);\r
+             N_VDestroy(cv_mem.cv_acorQ);\r
+             N_VDestroy(cv_mem.cv_yQ);\r
+             N_VDestroy(cv_mem.cv_tempvQ);\r
+             for (i=0; i < j; i++) N_VDestroy(cv_mem.cv_znQ[i]);\r
+             return(false);\r
+           }\r
+         }\r
+\r
+         /* Store the value of qmax used here */\r
+         cv_mem.cv_qmax_allocQ = cv_mem.cv_qmax;\r
+\r
+         /* Update solver workspace lengths */\r
+         cv_mem.cv_lrw += (cv_mem.cv_qmax + 5)*cv_mem.cv_lrw1Q;\r
+         cv_mem.cv_liw += (cv_mem.cv_qmax + 5)*cv_mem.cv_liw1Q;\r
+\r
+         return(true);\r
+       }\r
+\r
+       private int CVodeQuadSStolerances(CVodeMemRec cv_mem, double reltolQ, double abstolQ)\r
+       {\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES",\r
+                          "CVodeQuadSStolerances", MSGCV_NO_MEM);    \r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Ckeck if quadrature was initialized? */\r
+\r
+         if (cv_mem.cv_QuadMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_QUAD, "CVODES",\r
+                          "CVodeQuadSStolerances", MSGCV_NO_QUAD); \r
+           return(CV_NO_QUAD);\r
+         }\r
+\r
+         /* Test user-supplied tolerances */\r
+\r
+         if (reltolQ < ZERO) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeQuadSStolerances", MSGCV_BAD_RELTOLQ);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         if (abstolQ < 0) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
+                          "CVodeQuadSStolerances", MSGCV_BAD_ABSTOLQ);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Copy tolerances into memory */\r
+\r
+         cv_mem.cv_itolQ = CV_SS;\r
+\r
+         cv_mem.cv_reltolQ  = reltolQ;\r
+         cv_mem.cv_SabstolQ = abstolQ;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /*\r
+        * CVodeAdjInit\r
+        *\r
+        * This routine initializes ASA and allocates space for the adjoint \r
+        * memory structure.\r
+        */\r
+\r
+       private int CVodeAdjInit(CVodeMemRec cv_mem, int steps, int interp)\r
+       {\r
+         CVadjMemRec ca_mem;\r
+         int i, ii;\r
+\r
+         /* ---------------\r
+          * Check arguments\r
+          * --------------- */\r
+\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODEA", "CVodeAdjInit", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         if (steps <= 0) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeAdjInit", MSGCV_BAD_STEPS);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         if ( (interp != CV_HERMITE) && (interp != CV_POLYNOMIAL) ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeAdjInit", MSGCV_BAD_INTERP);\r
+           return(CV_ILL_INPUT);\r
+         } \r
+\r
+         /* ----------------------------\r
+          * Allocate CVODEA memory block\r
+          * ---------------------------- */\r
+\r
+         ca_mem = null;\r
+         ca_mem = new CVadjMemRec();\r
+\r
+         /* Attach ca_mem to CVodeMem structure */\r
+\r
+         cv_mem.cv_adj_mem = ca_mem;\r
+\r
+         /* ------------------------------\r
+          * Initialization of check points\r
+          * ------------------------------ */\r
+\r
+         /* Set Check Points linked list to NULL */\r
+         ca_mem.ck_mem = null;\r
+\r
+         /* Initialize nckpnts to ZERO */\r
+         ca_mem.ca_nckpnts = 0;\r
+\r
+         /* No interpolation data is available */\r
+         ca_mem.ca_ckpntData = null;\r
+\r
+         /* ------------------------------------\r
+          * Initialization of interpolation data\r
+          * ------------------------------------ */\r
+\r
+         /* Interpolation type */\r
+\r
+         ca_mem.ca_IMtype = interp;\r
+\r
+         /* Number of steps between check points */\r
+\r
+         ca_mem.ca_nsteps = steps;\r
+\r
+         /* Allocate space for the array of Data Point structures */\r
+\r
+         ca_mem.dt_mem = null;\r
+         ca_mem.dt_mem = new DtpntMemRec[steps+1];\r
+         if (ca_mem.dt_mem == null) {\r
+           ca_mem = null;\r
+           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);\r
+           return(CV_MEM_FAIL);\r
+         }\r
+\r
+         for (i=0; i<=steps; i++) { \r
+           ca_mem.dt_mem[i] = null;\r
+           ca_mem.dt_mem[i] = new DtpntMemRec();\r
+           if (ca_mem.dt_mem[i] == null) {\r
+             for(ii=0; ii<i; ii++) {ca_mem.dt_mem[ii] = null;}\r
+             ca_mem.dt_mem = null;\r
+             ca_mem = null;\r
+             cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);\r
+             return(CV_MEM_FAIL);\r
+           }\r
+         }\r
+\r
+         /* Attach functions for the appropriate interpolation module */\r
+         \r
+         switch(interp) {\r
+\r
+         case CV_HERMITE:\r
+           \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_select;\r
+\r
+           break;\r
+           \r
+         case CV_POLYNOMIAL:\r
+         \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_select;\r
+\r
+           break;\r
+\r
+         }\r
+\r
+         /* The interpolation module has not been initialized yet */\r
+\r
+         ca_mem.ca_IMmallocDone = false;\r
+\r
+         /* By default we will store but not interpolate sensitivities\r
+          *  - IMstoreSensi will be set in CVodeF to SUNFALSE if FSA is not enabled\r
+          *    or if the user can force this through CVodeSetAdjNoSensi \r
+          *  - IMinterpSensi will be set in CVodeB to SUNTRUE if IMstoreSensi is\r
+          *    SUNTRUE and if at least one backward problem requires sensitivities */\r
+\r
+         ca_mem.ca_IMstoreSensi = true;\r
+         ca_mem.ca_IMinterpSensi = false;\r
+\r
+         /* ------------------------------------\r
+          * Initialize list of backward problems\r
+          * ------------------------------------ */\r
+\r
+         ca_mem.cvB_mem = null;\r
+         ca_mem.ca_bckpbCrt = null;\r
+         ca_mem.ca_nbckpbs = 0;\r
+\r
+         /* --------------------------------\r
+          * CVodeF and CVodeB not called yet\r
+          * -------------------------------- */\r
+\r
+         ca_mem.ca_firstCVodeFcall = true;\r
+         ca_mem.ca_tstopCVodeFcall = false;\r
+\r
+         ca_mem.ca_firstCVodeBcall = true;\r
+\r
+         /* ---------------------------------------------\r
+          * ASA initialized and allocated\r
+          * --------------------------------------------- */\r
+\r
+         cv_mem.cv_adj = true;\r
+         cv_mem.cv_adjMallocDone = true;\r
+\r
+         return(CV_SUCCESS);\r
+       } \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
+               \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
+}
\ No newline at end of file