Started work on cvsRoberts_ASAi_dns.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 30 Mar 2018 22:16:06 +0000 (22:16 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 30 Mar 2018 22:16:06 +0000 (22:16 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15433 ba61647d-9d00-f842-95cd-605cb4296b96

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

index c1cc127..c8fb0c6 100644 (file)
@@ -2,6 +2,9 @@ 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
@@ -146,6 +149,13 @@ public abstract class CVODES {
        final int CV_CENTERED = 1;\r
        final int CV_FORWARD = 2;\r
        \r
+       \r
+        // CV_HERMITE specifies cubic Hermite interpolation.\r
+        // CV_POYNOMIAL specifies the polynomial interpolation\r
+       /* interp */\r
+       final int CV_HERMITE = 1;\r
+       final int CV_POLYNOMIAL = 2;\r
+\r
        /* \r
         * ----------------------------------------\r
         * CVODES return flags\r
@@ -575,7 +585,8 @@ public abstract class CVODES {
     final int cvsRoberts_dns_uw = 3;\r
     final int cvsRoberts_dnsL = 4;\r
     final int cvsRoberts_FSA_dns = 5;\r
-    int problem = cvsRoberts_dnsL;\r
+    final int cvsRoberts_ASAi_dns = 6;\r
+    int problem = cvsRoberts_FSA_dns;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -625,7 +636,19 @@ public abstract class CVODES {
          public int ewt(NVector y, NVector w, UserData user_data) {\r
               return 0;\r
           }\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
@@ -666,6 +689,9 @@ public abstract class CVODES {
            else if (problem == cvsRoberts_FSA_dns) {\r
                runcvsRoberts_FSA_dns();\r
            }\r
+           else if (problem == cvsRoberts_ASAi_dns) {\r
+               runcvsRoberts_ASAi_dns();\r
+           }\r
            else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                runcvsDirectDemo_Problem_1();\r
            }\r
@@ -1677,7 +1703,7 @@ public abstract class CVODES {
                        case 9:\r
                                System.out.println("Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995");\r
                                break;\r
-                       case 910:\r
+                       case 10:\r
                                System.out.println("Example manual at t = 4.0E8 y[0] = 5.2106E-6 y[1] = 2.0842E-11 y[2] = 0.99999");\r
                                break;\r
                        case 11:\r
@@ -1740,7 +1766,223 @@ 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
+               \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 allocting the internal \r
+        memory needed for backward integration.*/\r
+        steps = STEPS; /* no. of integration steps between two consecutive ckeckpoints*/\r
+        flag = CVodeAdjInit(cvode_mem, steps, CV_HERMITE);\r
+        /*\r
+        flag = CVodeAdjInit(cvode_mem, steps, CV_POLYNOMIAL);\r
+        */\r
+\r
+\r
+       }\r
 \r
        private void PrintFinalStats(CVodeMemRec cv_mem) {\r
            System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
@@ -2194,7 +2436,7 @@ public abstract class CVODES {
                    ydot[2] = 3.0E7*y[1]*y[1];\r
                    ydot[1] = -ydot[0] - ydot[2];\r
                }\r
-               else if (problem == cvsRoberts_FSA_dns) {\r
+               else if ((problem == cvsRoberts_FSA_dns) || (problem == cvsRoberts_ASAi_dns)) {\r
                        double p[] = user_data.array;\r
                ydot[0] = -p[0]*y[0] + p[1]*y[1]*y[2];\r
                ydot[2] = p[2]*y[1]*y[1];\r
@@ -2210,6 +2452,9 @@ public abstract class CVODES {
        public abstract int f(double t, NVector yv, NVector ydotv, UserData user_data);\r
        \r
        private int fQTestMode(double t, NVector x, NVector y, UserData user_data) {\r
+               if (problem == cvsRoberts_ASAi_dns) {\r
+                       y.data[0] = x.data[2];\r
+               }\r
                return 0;\r
        }\r
        \r
@@ -2293,7 +2538,7 @@ public abstract class CVODES {
                    J[2][1] = 6.0E7 * y[1];\r
                    J[2][2] = ZERO;\r
                }\r
-               else if (problem == cvsRoberts_FSA_dns) {\r
+               else if ((problem == cvsRoberts_FSA_dns) || (problem == cvsRoberts_ASAi_dns)) {\r
                        double p[] = data.array;\r
                        J[0][0] = -p[0];\r
                        J[0][1] = p[1]*y[2];\r
@@ -2320,7 +2565,7 @@ public abstract class CVODES {
        public abstract int Jac(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, NVector tmp2, NVector tmp3);\r
        \r
        public int ewtTestMode(NVector y, NVector w, UserData user_data) {\r
-               if ((problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_dns)) {\r
+               if ((problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_dns) || (problem == cvsRoberts_ASAi_dns)) {\r
                        //final double RTOL = 1.0E-4; // scalar relative tolerance\r
                        final double RTOL = 1.0E-12;\r
                        //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
@@ -2352,7 +2597,57 @@ public abstract class CVODES {
        \r
        public abstract int ewt(NVector y, NVector w, UserData user_data);\r
        \r
+       public int fBTestMode(double t, NVector y, NVector yB, NVector yBdot, UserData user_dataB) {\r
+               if (problem == cvsRoberts_ASAi_dns) {\r
+                       double p[] = user_dataB.array;\r
+                       double l10,l21,y12;\r
+                       l10 = yB.data[1] - yB.data[0];\r
+                       l21 = yB.data[2] - yB.data[1];\r
+                       y12 = y.data[1] * y.data[2];\r
+                       yBdot.data[0] = -p[0] * l10;\r
+                       yBdot.data[1] = p[1]*y.data[2]*l10 - 2.0*p[2]*y.data[1]*l21;\r
+                       yBdot.data[2] = p[1]*y.data[1]*l10 - 1.0;\r
+               }\r
+               return 0;\r
+       }\r
+       \r
+       public abstract int fB(double t, NVector y, NVector yB, NVector yBdot, UserData user_dataB);\r
+       \r
+       public int JacBTestMode(double t, NVector y, NVector yB, NVector fyB, double JB[][], UserData user_dataB,\r
+                       NVector tmp1B, NVector tmp2B, NVector tmp3B) {\r
+               if (problem == cvsRoberts_ASAi_dns) {\r
+                       double p[] = user_dataB.array;\r
+                       JB[0][0] = p[0];\r
+                       JB[0][1] = -p[0];\r
+                       JB[0][2] = ZERO;\r
+                       JB[1][0] = -p[1]*y.data[2];\r
+                       JB[1][1] = p[1]*y.data[2] + 2.0*p[2]*y.data[1];\r
+                       JB[1][2] = -2.0*p[2]*y.data[1];\r
+                       JB[2][0] = -p[1]*y.data[1];\r
+                       JB[2][1] = p[1]*y.data[1];\r
+                       JB[2][2] = ZERO;\r
+               }\r
+               return 0;\r
+       }\r
        \r
+       public abstract int JacB(double t, NVector y, NVector yB, NVector fyB, double JB[][], UserData user_dataB,\r
+                       NVector tmp1B, NVector tmp2B, NVector tmp3B);\r
+       \r
+       // fQB routine.  Computes integrand for quadratures\r
+       public int fQBTestMode(double t, NVector y, NVector yB, NVector qBdot, UserData user_dataB) {\r
+               if (problem == cvsRoberts_ASAi_dns) {\r
+                       double l10, l21, y12;\r
+                       l10 = yB.data[1] - yB.data[0];\r
+                       l21 = yB.data[2] - yB.data[1];\r
+                       y12 = y.data[1] * y.data[2];\r
+                       qBdot.data[0] = y.data[0] * l10;\r
+                       qBdot.data[1] = -y12 * l10;\r
+                       qBdot.data[2] = y.data[1] *y.data[1] * l21;\r
+               }\r
+               return 0;\r
+       }\r
+       \r
+       public abstract int fQB(double t, NVector y, NVector yB, NVector qBdot, UserData user_dataB);\r
        \r
        // Types: struct CVodeMemRec, CVodeMem\r
        // -----------------------------------------------------------------\r
@@ -2410,6 +2705,7 @@ public abstract class CVODES {
          boolean cv_quadr;       /* SUNTRUE if integrating quadratures            */\r
 \r
          //CVQuadRhsFn cv_fQ;          /* q' = fQ(t, y(t))                              */\r
+         int cv_fQ;\r
 \r
          boolean cv_errconQ;     /* SUNTRUE if quadrs. are included in error test */\r
 \r
@@ -2634,11 +2930,11 @@ public abstract class CVODES {
            Space requirements for CVODES \r
            -----------------------------*/\r
 \r
-         long cv_lrw1;        /* no. of realtype words in 1 N_Vector y           */ \r
+         long cv_lrw1;        /* no. of double words in 1 N_Vector y           */ \r
          long cv_liw1;        /* no. of integer words in 1 N_Vector y            */ \r
-         long cv_lrw1Q;       /* no. of realtype words in 1 N_Vector yQ          */ \r
+         long cv_lrw1Q;       /* no. of double words in 1 N_Vector yQ          */ \r
          long cv_liw1Q;       /* no. of integer words in 1 N_Vector yQ           */ \r
-         long cv_lrw;             /* no. of realtype words in CVODES work vectors    */\r
+         long cv_lrw;             /* no. of double words in CVODES work vectors    */\r
          long cv_liw;             /* no. of integer words in CVODES work vectors     */\r
 \r
          /*----------------\r
@@ -2762,7 +3058,7 @@ public abstract class CVODES {
 \r
          boolean cv_adj;             /* SUNTRUE if performing ASA                */\r
 \r
-         CVodeMemRec cv_adj_mem; /* Pointer to adjoint memory structure      */\r
+         CVadjMemRec cv_adj_mem; /* Pointer to adjoint memory structure      */\r
 \r
          boolean cv_adjMallocDone;\r
 \r
@@ -11707,6 +12003,560 @@ 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
+\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
+         if (ca_mem == null) {\r
+           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);\r
+           return(CV_MEM_FAIL);\r
+         }\r
+\r
+         /* Attach ca_mem to CVodeMem structure */\r
+\r
+         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;\r
+           //ca_mem.ca_IMfree   = CVAhermiteFree;\r
+           //ca_mem.ca_IMget    = CVAhermiteGetY;\r
+           //ca_mem.ca_IMstore  = CVAhermiteStorePnt;\r
+\r
+           break;\r
+           \r
+         case CV_POLYNOMIAL:\r
+         \r
+           //ca_mem.ca_IMmalloc = CVApolynomialMalloc;\r
+           //ca_mem.ca_IMfree   = CVApolynomialFree;\r
+           //ca_mem.ca_IMget    = CVApolynomialGetY;\r
+           //ca_mem.ca_IMstore  = CVApolynomialStorePnt;\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
+                 //cvaIMMallocFn   ca_IMmalloc; \r
+                 //cvaIMFreeFn     ca_IMfree;\r
+                 //cvaIMStorePntFn 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
+\r
+                 /* Right hand side function for backward run */\r
+                 //CVRhsFnB cv_f;\r
+                 //CVRhsFnBS cv_fs;\r
+\r
+                 /* Right hand side quadrature function for backward run */\r
+                 //CVQuadRhsFnB cv_fQ;\r
+                 //CVQuadRhsFnBS cv_fQs;\r
+\r
+                 /* User user_data */\r
+                 UserData cv_user_data;\r
+                   \r
+                 /* Memory block for a linear solver's interface to CVODEA */\r
+                 //void *cv_lmem;\r
+\r
+                 /* Function to free any memory allocated by the linear solver */\r
+                 //int (*cv_lfree)(CVodeBMem cvB_mem);\r
+\r
+                 /* Memory block for a preconditioner's module interface to CVODEA */ \r
+                 //void *cv_pmem;\r
+\r
+                 /* Function to free any memory allocated by the preconditioner module */\r
+                 //int (*cv_pfree)(CVodeBMem cvB_mem);\r
+\r
+                 /* Time at which to extract solution / quadratures */\r
+                 double cv_tout;\r
+                 \r
+                 /* Workspace Nvector */\r
+                 NVector cv_y;\r
+\r
+                 /* Pointer to next structure in list */\r
+                 CVodeBMemRec cv_next;\r
+\r
+               };\r
+               \r
+               /*\r
+                * -----------------------------------------------------------------\r
+                * Types : struct CkpntMemRec, CkpntMem\r
+                * -----------------------------------------------------------------\r
+                * The type CkpntMem is type pointer to struct CkpntMemRec.\r
+                * This structure contains fields to store all information at a\r
+                * check point that is needed to 'hot' start cvodes.\r
+                * -----------------------------------------------------------------\r
+                */\r
+\r
+               class CkpntMemRec {\r
+\r
+                 /* Integration limits */\r
+                 double ck_t0;\r
+                 double ck_t1;\r
+                   \r
+                 /* Nordsieck History Array */\r
+                 NVector ck_zn[] = new NVector[L_MAX];\r
+                   \r
+                 /* Do we need to carry quadratures? */\r
+                 boolean ck_quadr;\r
+                   \r
+                 /* Nordsieck History Array for quadratures */\r
+                 NVector ck_znQ[] = new NVector[L_MAX];\r
+                   \r
+                 /* Do we need to carry sensitivities? */\r
+                 boolean ck_sensi;\r
+\r
+                 /* number of sensitivities */\r
+                 int ck_Ns;\r
+\r
+                 /* Nordsieck History Array for sensitivities */\r
+                 NVector ck_znS[][] = new NVector[L_MAX][];\r
+\r
+                 /* Do we need to carry quadrature sensitivities? */\r
+                 boolean ck_quadr_sensi;\r
+\r
+                 /* Nordsieck History Array for quadrature sensitivities */\r
+                 NVector ck_znQS[][] = new NVector[L_MAX][];\r
+                   \r
+                 /* Was ck_zn[qmax] allocated?\r
+                    ck_zqm = 0    - no\r
+                    ck_zqm = qmax - yes      */\r
+                 int ck_zqm;\r
+                   \r
+                 /* Step data */\r
+                 long ck_nst;\r
+                 double ck_tretlast;\r
+                 int      ck_q;\r
+                 int      ck_qprime;\r
+                 int      ck_qwait;\r
+                 int      ck_L;\r
+                 double ck_gammap;\r
+                 double ck_h;\r
+                 double ck_hprime;\r
+                 double ck_hscale;\r
+                 double ck_eta;\r
+                 double ck_etamax;\r
+                 double ck_tau[] = new double[L_MAX+1];\r
+                 double ck_tq[] = new double[NUM_TESTS+1];\r
+                 double ck_l[] = new double[L_MAX];\r
+                   \r
+                 /* Saved values */\r
+                 double ck_saved_tq5;\r
+                   \r
+                 /* Pointer t next structure in list */\r
+                 CkpntMemRec ck_next;\r
+                   \r
+               };\r
+               \r
+               /*\r
+                * -----------------------------------------------------------------\r
+                * Type : struct DtpntMemRec\r
+                * -----------------------------------------------------------------\r
+                * This structure contains fields to store all information at a\r
+                * data point that is needed to interpolate solution of forward\r
+                * simulations. Its content field depends on IMtype.\r
+                * -----------------------------------------------------------------\r
+                */\r
+                 \r
+               class DtpntMemRec {\r
+                 double t;    /* time */\r
+                 //void *content; /* IMtype-dependent content */\r
+               };\r
 \r
 \r
 \r