CVODES split int CVODES and CVODES_ASA for running runcvsRoberts_ASAi_dns().
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 18:24:33 +0000 (18:24 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 3 Apr 2018 18:24:33 +0000 (18:24 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15435 ba61647d-9d00-f842-95cd-605cb4296b96

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

index da53192..756e15f 100644 (file)
@@ -2,9 +2,6 @@ package gov.nih.mipav.model.algorithms;
 \r
 import java.io.RandomAccessFile;\r
 \r
-import gov.nih.mipav.model.algorithms.CVODES.CVodeMemRec;\r
-import gov.nih.mipav.model.algorithms.CVODES.NVector;\r
-import gov.nih.mipav.model.algorithms.CVODES.UserData;\r
 import gov.nih.mipav.model.structures.jama.LinearEquations2;\r
 import gov.nih.mipav.view.MipavUtil;\r
 \r
@@ -605,7 +602,7 @@ public abstract class CVODES {
     final int cvsRoberts_dnsL = 4;\r
     final int cvsRoberts_FSA_dns = 5;\r
     final int cvsRoberts_ASAi_dns = 6;\r
-    int problem = cvsRoberts_FSA_dns;\r
+    int problem = cvsRoberts_ASAi_dns;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -624,9 +621,11 @@ public abstract class CVODES {
        final double P1_DTOUT = 2.214773875;\r
        final double P1_TOL_FACTOR = 1.0E4;\r
        \r
-       // USe the following code to call CVODES from another module:\r
+       // Use the following code to call CVODES or CVODES_ASA from another module:\r
        /*boolean testme = true;\r
-    class CVODEStest extends CVODES {\r
+       Use only one of the following 2 lines\r
+    class CVODEStest extends CVODES { // if not running runcvsRoberts_ASAi_dns()\r
+    class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns()\r
          public CVODEStest() {\r
                  super();\r
          }\r
@@ -710,7 +709,7 @@ public abstract class CVODES {
            }\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
+               // Not run here.  Run in CVODES_ASA.\r
            //}\r
            else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                runcvsDirectDemo_Problem_1();\r
@@ -1786,237 +1785,7 @@ public abstract class CVODES {
                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
        private void PrintFinalStats(CVodeMemRec cv_mem) {\r
            System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
@@ -2447,7 +2216,7 @@ public abstract class CVODES {
            System.out.println("Error overrun = " + ero);\r
        }\r
        \r
-       private void N_VNew_Serial(NVector y, int length) {\r
+       protected void N_VNew_Serial(NVector y, int length) {\r
            if ((y != null) && (length > 0)) {\r
                y.data = new double[length];\r
                y.own_data = true;\r
@@ -2462,7 +2231,7 @@ public abstract class CVODES {
         * @param user_data\r
         * @return\r
         */\r
-       private int fTestMode(double t, NVector yv, NVector ydotv, UserData user_data) {\r
+       protected int fTestMode(double t, NVector yv, NVector ydotv, UserData user_data) {\r
                double y[] = yv.getData();\r
                double ydot[] = ydotv.getData();\r
                if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_dnsL)) {\r
@@ -3120,7 +2889,7 @@ public abstract class CVODES {
         // If an initialization error occurs, CVodeCreate prints an error\r
         // message to standard err and returns NULL.\r
 \r
-     private CVodeMemRec CVodeCreate(int lmm, int iter) {\r
+     protected CVodeMemRec CVodeCreate(int lmm, int iter) {\r
         int maxord;\r
          CVodeMemRec cv_mem;\r
 \r
@@ -3554,14 +3323,14 @@ public abstract class CVODES {
        return true;\r
      }\r
      \r
-     private NVector N_VClone(NVector y) {\r
+     protected NVector N_VClone(NVector y) {\r
          NVector x = new NVector();\r
          x.data = y.data.clone();\r
          x.own_data = y.own_data;\r
          return x;\r
      }\r
 \r
-     private void N_VDestroy(NVector x) {\r
+     protected void N_VDestroy(NVector x) {\r
         x.data = null;\r
         x = null;\r
      }\r
@@ -3617,7 +3386,7 @@ public abstract class CVODES {
        return(CV_SUCCESS);\r
      }\r
      \r
-     private int CVodeWFtolerances(CVodeMemRec cv_mem, int efun_select)\r
+     protected int CVodeWFtolerances(CVodeMemRec cv_mem, int efun_select)\r
        {\r
 \r
          if (cv_mem== null) {\r
@@ -3828,7 +3597,7 @@ public abstract class CVODES {
          SUNLinearSolver_Type type;\r
      }\r
      \r
-     private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
+     protected SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
         int matrixRows = A.length;\r
         int matrixColumns = A[0].length;\r
         if (matrixRows != matrixColumns) {\r
@@ -5071,7 +4840,7 @@ public abstract class CVODES {
   }\r
 \r
 \r
-  private void N_VScale_Serial(double c, NVector x, NVector z)\r
+  protected void N_VScale_Serial(double c, NVector x, NVector z)\r
   {\r
     int i, N;\r
     double xd[], zd[];\r
@@ -6483,7 +6252,7 @@ else                return(snrm);
   * may also be called directly by the user.\r
   */\r
 \r
- private int CVodeGetDky(CVodeMemRec cv_mem, double t, int k, NVector dky)\r
+ protected int CVodeGetDky(CVodeMemRec cv_mem, double t, int k, NVector dky)\r
  {\r
    double s, c, r;\r
    double tfuzz, tp, tn1;\r
@@ -11140,7 +10909,7 @@ else                return(snrm);
            * Specifies the maximum number of integration steps\r
            */\r
 \r
-          private int CVodeSetMaxNumSteps(CVodeMemRec cv_mem, long mxsteps)\r
+          protected int CVodeSetMaxNumSteps(CVodeMemRec cv_mem, long mxsteps)\r
           {\r
 \r
             if (cv_mem==null) {\r
@@ -11164,7 +10933,7 @@ else                return(snrm);
            * step try.\r
            */\r
 \r
-          private int CVodeSetMaxErrTestFails(CVodeMemRec cv_mem, int maxnef)\r
+          protected int CVodeSetMaxErrTestFails(CVodeMemRec cv_mem, int maxnef)\r
           {\r
 \r
             if (cv_mem==null) {\r
@@ -11489,7 +11258,7 @@ else                return(snrm);
         * Function to create an array of new serial vectors. \r
         */\r
 \r
-       private NVector[] N_VCloneVectorArray_Serial(int count, NVector w)\r
+       protected NVector[] N_VCloneVectorArray_Serial(int count, NVector w)\r
        {\r
          NVector vs[];\r
          int j;\r
@@ -11511,7 +11280,7 @@ else                return(snrm);
          return(vs);\r
        }\r
        \r
-       private NVector N_VClone_Serial(NVector w)\r
+       protected NVector N_VClone_Serial(NVector w)\r
        {\r
          NVector v;\r
          double data[];\r
@@ -11543,7 +11312,7 @@ else                return(snrm);
         * Function to free an array created with N_VCloneVectorArray_Serial\r
         */\r
 \r
-       private void N_VDestroyVectorArray_Serial(NVector vs[], int count)\r
+       protected void N_VDestroyVectorArray_Serial(NVector vs[], int count)\r
        {\r
          int j;\r
 \r
@@ -12037,437 +11806,119 @@ else                return(snrm);
          \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
+       class CVadjMemRec {\r
+           \r
+                 /* --------------------\r
+                  * Forward problem data\r
+                  * -------------------- */\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
+                 /* Integration interval */\r
+                 double ca_tinitial, ca_tfinal;\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
+                 /* Flag for first call to CVodeF */\r
+                 boolean ca_firstCVodeFcall;\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
+                 /* Flag if CVodeF was called with TSTOP */\r
+                 boolean ca_tstopCVodeFcall;\r
+                 double ca_tstopCVodeF;\r
+                   \r
+                 /* ----------------------\r
+                  * Backward problems data\r
+                  * ---------------------- */\r
 \r
-         /* Initialize znQ[0] in the history array */\r
-         N_VScale_Serial(ONE, yQ0, cv_mem.cv_znQ[0]);\r
+                 /* Storage for backward problems */\r
+                 CVodeBMemRec cvB_mem;\r
 \r
-         /* Copy the input parameters into CVODES state */\r
-         cv_mem.cv_fQ = fQ;\r
+                 /* Number of backward problems */\r
+                 int ca_nbckpbs;\r
 \r
-         /* Initialize counters */\r
-         cv_mem.cv_nfQe  = 0;\r
-         cv_mem.cv_netfQ[0] = 0;\r
+                 /* Address of current backward problem */\r
+                 CVodeBMemRec ca_bckpbCrt;\r
 \r
-         /* Quadrature integration turned ON */\r
-         cv_mem.cv_quadr = true;\r
-         cv_mem.cv_QuadMallocDone = true;\r
+                 /* Flag for first call to CVodeB */\r
+                 boolean ca_firstCVodeBcall;\r
+                   \r
+                 /* ----------------\r
+                  * Check point data\r
+                  * ---------------- */\r
 \r
-         /* Quadrature initialization was successfull */\r
-         return(CV_SUCCESS);\r
-       }\r
+                 /* Storage for check point information */\r
+                 CkpntMemRec ck_mem;\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
+                 /* Number of check points */\r
+                 int ca_nckpnts;\r
 \r
-       private boolean cvQuadAllocVectors(CVodeMemRec cv_mem, NVector tmpl) \r
-       {\r
-         int i, j;\r
+                 /* address of the check point structure for which data is available */\r
+                 CkpntMemRec ca_ckpntData;\r
+                   \r
+                 /* ------------------\r
+                  * Interpolation data\r
+                  * ------------------ */\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
+                 /* Number of steps between 2 check points */\r
+                 long ca_nsteps;\r
+                 \r
+                 /* Storage for data from forward runs */\r
+                 DtpntMemRec dt_mem[];\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
+                 /* Actual number of data points in dt_mem (typically np=nsteps+1) */\r
+                 long ca_np;\r
+                   \r
+                 /* Interpolation type */\r
+                 int ca_IMtype;\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
+                 /* Functions set by the interpolation module */\r
+                 int   ca_IMmalloc; \r
+                 //cvaIMFreeFn     ca_IMfree;\r
+                 int ca_IMstore; /* store a new interpolation point */\r
+                 //cvaIMGetYFn     ca_IMget;   /* interpolate forward solution    */\r
 \r
-         /* Allocate zQn[0] ... zQn[maxord] */\r
+                 /* Flags controlling the interpolation module */\r
+                 boolean ca_IMmallocDone;   /* IM initialized? */\r
+                 boolean ca_IMnewData;      /* new data available in dt_mem?*/\r
+                 boolean ca_IMstoreSensi;   /* store sensitivities? */\r
+                 boolean ca_IMinterpSensi;  /* interpolate sensitivities? */\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
+                 /* Workspace for the interpolation module */\r
+                 NVector ca_Y[] = new NVector[L_MAX];     /* pointers to zn[i] */\r
+                 NVector ca_YS[][] = new NVector[L_MAX][];   /* pointers to znS[i] */\r
+                 double ca_T[] = new double[L_MAX];\r
 \r
-         /* Store the value of qmax used here */\r
-         cv_mem.cv_qmax_allocQ = cv_mem.cv_qmax;\r
+                 /* -------------------------------\r
+                  * Workspace for wrapper functions\r
+                  * ------------------------------- */\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
+                 NVector ca_ytmp;\r
 \r
-         return(true);\r
-       }\r
+                 NVector ca_yStmp[];\r
+                   \r
+               };\r
 \r
-       private int CVodeQuadSStolerances(CVodeMemRec cv_mem, double reltolQ, double abstolQ)\r
-       {\r
+               /*\r
+                * -----------------------------------------------------------------\r
+                * Type : struct CVodeBMemRec\r
+                * -----------------------------------------------------------------\r
+                * The type CVodeBMem is a pointer to a structure which stores all\r
+                * information for ONE backward problem.\r
+                * The CVadjMem structure contains a linked list of CVodeBMem pointers\r
+                * -----------------------------------------------------------------\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
+               class CVodeBMemRec {\r
 \r
-         /* Ckeck if quadrature was initialized? */\r
+                 /* Index of this backward problem */\r
+                 int cv_index;\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
+                 /* Time at which the backward problem is initialized */\r
+                 double cv_t0;\r
+                 \r
+                 /* CVODES memory for this backward problem */\r
+                 CVodeMemRec cv_mem;\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
-       class CVadjMemRec {\r
-           \r
-                 /* --------------------\r
-                  * Forward problem data\r
-                  * -------------------- */\r
-\r
-                 /* Integration interval */\r
-                 double ca_tinitial, ca_tfinal;\r
-\r
-                 /* Flag for first call to CVodeF */\r
-                 boolean ca_firstCVodeFcall;\r
-\r
-                 /* Flag if CVodeF was called with TSTOP */\r
-                 boolean ca_tstopCVodeFcall;\r
-                 double ca_tstopCVodeF;\r
-                   \r
-                 /* ----------------------\r
-                  * Backward problems data\r
-                  * ---------------------- */\r
-\r
-                 /* Storage for backward problems */\r
-                 CVodeBMemRec cvB_mem;\r
-\r
-                 /* Number of backward problems */\r
-                 int ca_nbckpbs;\r
-\r
-                 /* Address of current backward problem */\r
-                 CVodeBMemRec ca_bckpbCrt;\r
-\r
-                 /* Flag for first call to CVodeB */\r
-                 boolean ca_firstCVodeBcall;\r
-                   \r
-                 /* ----------------\r
-                  * Check point data\r
-                  * ---------------- */\r
-\r
-                 /* Storage for check point information */\r
-                 CkpntMemRec ck_mem;\r
-\r
-                 /* Number of check points */\r
-                 int ca_nckpnts;\r
-\r
-                 /* address of the check point structure for which data is available */\r
-                 CkpntMemRec ca_ckpntData;\r
-                   \r
-                 /* ------------------\r
-                  * Interpolation data\r
-                  * ------------------ */\r
-\r
-                 /* Number of steps between 2 check points */\r
-                 long ca_nsteps;\r
-                 \r
-                 /* Storage for data from forward runs */\r
-                 DtpntMemRec dt_mem[];\r
-\r
-                 /* Actual number of data points in dt_mem (typically np=nsteps+1) */\r
-                 long ca_np;\r
-                   \r
-                 /* Interpolation type */\r
-                 int ca_IMtype;\r
-\r
-                 /* Functions set by the interpolation module */\r
-                 int   ca_IMmalloc; \r
-                 //cvaIMFreeFn     ca_IMfree;\r
-                 int ca_IMstore; /* store a new interpolation point */\r
-                 //cvaIMGetYFn     ca_IMget;   /* interpolate forward solution    */\r
-\r
-                 /* Flags controlling the interpolation module */\r
-                 boolean ca_IMmallocDone;   /* IM initialized? */\r
-                 boolean ca_IMnewData;      /* new data available in dt_mem?*/\r
-                 boolean ca_IMstoreSensi;   /* store sensitivities? */\r
-                 boolean ca_IMinterpSensi;  /* interpolate sensitivities? */\r
-\r
-                 /* Workspace for the interpolation module */\r
-                 NVector ca_Y[] = new NVector[L_MAX];     /* pointers to zn[i] */\r
-                 NVector ca_YS[][] = new NVector[L_MAX][];   /* pointers to znS[i] */\r
-                 double ca_T[] = new double[L_MAX];\r
-\r
-                 /* -------------------------------\r
-                  * Workspace for wrapper functions\r
-                  * ------------------------------- */\r
-\r
-                 NVector ca_ytmp;\r
-\r
-                 NVector ca_yStmp[];\r
-                   \r
-               };\r
-\r
-               /*\r
-                * -----------------------------------------------------------------\r
-                * Type : struct CVodeBMemRec\r
-                * -----------------------------------------------------------------\r
-                * The type CVodeBMem is a pointer to a structure which stores all\r
-                * information for ONE backward problem.\r
-                * The CVadjMem structure contains a linked list of CVodeBMem pointers\r
-                * -----------------------------------------------------------------\r
-                */\r
-\r
-               class CVodeBMemRec {\r
-\r
-                 /* Index of this backward problem */\r
-                 int cv_index;\r
-\r
-                 /* Time at which the backward problem is initialized */\r
-                 double cv_t0;\r
-                 \r
-                 /* CVODES memory for this backward problem */\r
-                 CVodeMemRec cv_mem;\r
-\r
-                 /* Flags to indicate that this backward problem's RHS or quad RHS\r
-                  * require forward sensitivities */\r
-                 boolean cv_f_withSensi;\r
-                 boolean cv_fQ_withSensi;\r
+                 /* Flags to indicate that this backward problem's RHS or quad RHS\r
+                  * require forward sensitivities */\r
+                 boolean cv_f_withSensi;\r
+                 boolean cv_fQ_withSensi;\r
 \r
                  /* Right hand side function for backward run */\r
                  //CVRhsFnB cv_f;\r
@@ -12589,445 +12040,7 @@ else                return(snrm);
                  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
@@ -13042,391 +12055,6 @@ else                return(snrm);
                  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