Working on cvsRoberts_FSA_dns.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 14 Mar 2018 21:52:41 +0000 (21:52 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 14 Mar 2018 21:52:41 +0000 (21:52 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15416 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 8f40ef0..7e0392c 100644 (file)
@@ -563,9 +563,10 @@ public abstract class CVODES {
     final int cvDlsSetup_select = 1;\r
 \r
     final int cvsRoberts_dns = 1;\r
-    final int cvSDirectDemo_ls_Problem_1 = 2;\r
+    final int cvsDirectDemo_ls_Problem_1 = 2;\r
     final int cvsRoberts_dns_uw = 3;\r
-    int problem = cvsRoberts_dns_uw;\r
+    final int cvsRoberts_FSA_dns = 4;\r
+    int problem = cvsDirectDemo_ls_Problem_1;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -590,24 +591,24 @@ public abstract class CVODES {
          public CVODEStest() {\r
                  super();\r
          }\r
-         public int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
+         public int f(double t, NVector yv, NVector ydotv, UserData user_data) {\r
                  return 0;\r
          }\r
          \r
-         public int g(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
+         public int g(double t, NVector yv, double gout[], UserData user_data) {\r
                  return 0;\r
          }\r
          \r
-         public int fQ(double t, NVector x, NVector y, CVodeMemRec user_data) {\r
+         public int fQ(double t, NVector x, NVector y, UserData user_data) {\r
                  return 0;\r
          }\r
          \r
-         public int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, \r
+         public int Jac(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, \r
                          NVector tmp2, NVector tmp3) {\r
                  return 0;\r
          }\r
          \r
-         public int ewt(NVector y, NVector w, CVodeMemRec user_data) {\r
+         public int ewt(NVector y, NVector w, UserData user_data) {\r
               return 0;\r
           }\r
        \r
@@ -645,7 +646,10 @@ public abstract class CVODES {
            else if (problem == cvsRoberts_dns_uw) {\r
                runcvsRoberts_dns_uw();\r
            }\r
-           else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+           else if (problem == cvsRoberts_FSA_dns) {\r
+               runcvsRoberts_FSA_dns();\r
+           }\r
+           else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                runcvsDirectDemo_Problem_1();\r
            }\r
            return;\r
@@ -1156,6 +1160,280 @@ public abstract class CVODES {
                return;\r
        }\r
        \r
+       private void runcvsRoberts_FSA_dns() {\r
+               /* -----------------------------------------------------------------\r
+                * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, and\r
+                *                Radu Serban @ LLNL\r
+                * -----------------------------------------------------------------\r
+                * Example problem:\r
+                *\r
+                * The following is a simple example problem, with the coding\r
+                * needed for its solution by CVODES for Forward Sensitivity \r
+                * Analysis. The problem is from chemical kinetics, and consists\r
+                * 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: p1=0.04,\r
+                * p2=1e4, and p3=3e7. The problem is stiff.\r
+                * This program solves the problem with the BDF method, Newton\r
+                * iteration with the CVODES dense linear solver, and a\r
+                * user-supplied 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 the\r
+                * problem parameters p1, p2, and p3.\r
+                * The sensitivity right hand side is given analytically through the\r
+                * user routine fS (of type SensRhs1Fn).\r
+                * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and\r
+                * STAGGERED1) can be used and sensitivities may be included in the\r
+                * error test or not (error control set on SUNTRUE or SUNFALSE,\r
+                * respectively).\r
+                *\r
+                * Execution:\r
+                *\r
+                * If no sensitivities are desired:\r
+                *    % cvsRoberts_FSA_dns -nosensi\r
+                * If sensitivities are to be computed:\r
+                *    % cvsRoberts_FSA_dns -sensi sensi_meth err_con\r
+                * where sensi_meth is one of {sim, stg, stg1} and err_con is one of\r
+                * {t, f}.\r
+                * -----------------------------------------------------------------*/\r
+\r
+               /** Problem Constants */\r
+               final int NEQ = 3; // Number of equations\r
+               final int NS = 3; // Number of sensitivities computed\r
+               final double Y1 = 1.0; // Initial y components\r
+               final double Y2 = 0.0;\r
+               final double Y3 = 0.0;\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
+               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 T0 = 0.0; // initial time\r
+               final double T1 = 0.4; // first output time\r
+               final double TMULT = 10.0; // output time factor\r
+               //final int NOUT = 12; // number of output times\r
+               final int NOUT = 12;\r
+               int f = cvsRoberts_dns;\r
+               int g = cvsRoberts_dns;\r
+               int Jac = cvsRoberts_dns;\r
+               int ewt_select = cvEwtUser_select1;\r
+               // User data structure\r
+               double p[] = new double[]{0.04,1.0E4,3.0E7};\r
+               UserData data = new UserData();\r
+               data.array = p;\r
+               \r
+               // Set initial conditions\r
+               NVector y = new NVector();\r
+               NVector yS[];\r
+               NVector abstol = new NVector();\r
+               double yr[] = new double[]{Y1,Y2,Y3};\r
+               N_VNew_Serial(y, NEQ);\r
+               y.setData(yr);\r
+               double reltol = RTOL; // Set the scalar relative tolerance\r
+               // Set the vector absolute tolerance\r
+               double abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               N_VNew_Serial(abstol,NEQ);\r
+               abstol.setData(abstolr);\r
+               CVodeMemRec cvode_mem;\r
+               int flag;\r
+               int flagr;\r
+               double A[][];\r
+               SUNLinearSolver LS;\r
+               double tout;\r
+               int iout;\r
+               double t[] = new double[1];\r
+               int rootsfound[] = new int[2];\r
+               int i;\r
+               boolean sensi = true;\r
+               double pbar[] = new double[3];\r
+               \r
+               // Call CVodeCreate to create the solver memory and specify the\r
+               // Backward Differentiation Formula and the use of a Newton\r
+               // 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
+               // In loop, call CVode, print results, and test for error.\r
+               // Break out of loop when NOUT preset output times have been reached.\r
+               System.out.println(" \n3-species kinetics problem\n");\r
+               \r
+               // Sensitivity-related settings\r
+               if (sensi) {\r
+                       \r
+                       // Set parameter scaling factor\r
+                       pbar[0] = data.array[0];\r
+                       pbar[1] = data.array[1];\r
+                       pbar[2] = data.array[2];\r
+                       \r
+                       // Set sensitivity initial conditions\r
+                       yS = N_VCloneVectorArray_Serial(NS, y);\r
+                       if (yS == null) {\r
+                               return;\r
+                       }\r
+               } // if (sensi)\r
+               \r
+               iout = 0;\r
+               tout = T1;\r
+               \r
+               while (true) {\r
+                       flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       switch (iout) {\r
+                       case 0:\r
+                               System.out.println("Example manual: at t = 0.4 y[0] = 0.98517 y[1] = 3.3864E-5 y[2] = 1.4794E-2");\r
+                               break;\r
+                       case 1:\r
+                               System.out.println("Example manual: at t = 4.0 y[0] = 0.90552 y[1] = 2.2405E-5 y[2] = 9.4459E-2");\r
+                               break;\r
+                       case 2:\r
+                               System.out.println("Example manual: at t = 40.0 y[0] = 0.71583 y[1] = 9.1856E-6 y[2] = 0.28416");\r
+                               break;\r
+                       case 3:\r
+                               System.out.println("Example manual: at t = 400.0 y[0] = 0.45052 y[1] = 3.2229E-6 y[2] = 0.54947");\r
+                               break;\r
+                       case 4:\r
+                               System.out.println("Example manual: at t = 4.0E3 y[0] = 0.18317 y[1] = 8.9403E-7 y[2] = 0.81683");\r
+                               break;\r
+                       case 5:\r
+                               System.out.println("Example manual: at t = 4.0E4 y[0] = 3.8977E-2 y[1] = 1.6215E-7 y[2] = 0.96102");\r
+                               break;\r
+                       case 6:\r
+                               System.out.println("Example manual: at t = 4.0E5 y[0] = 4.9387E-3 y[1] = 1.9852E-8 y[2] = 0.99506");\r
+                               break;\r
+                       case 7:\r
+                               System.out.println("Example manual: at t = 4.0E6 y[0] = 5.1684E-4 y[1] = 2.0684E-9 y[2] = 0.99948");\r
+                               break;\r
+                       case 8:\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 9:\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 10:\r
+                               System.out.println("Example manual at t = 4.0E9 y[0] = 5.1881E-7 y[1] = 2.0752E-12 y[2] = 1.0000");\r
+                               break;\r
+                       case 11:\r
+                               System.out.println("Example manual at t = 4.0E10 y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000");\r
+                               break;\r
+                       }\r
+                       System.out.println("MIPAV: at t = " + t[0] + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
+                       \r
+                       if (flag == CV_ROOT_RETURN) {\r
+                           flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
+                           if (flagr == CV_MEM_NULL) {\r
+                               return;\r
+                           }\r
+                           System.out.println("Roots found are " + rootsfound[0] + " and " + rootsfound[1]);\r
+                       }\r
+                       \r
+                       if (flag < 0) {\r
+                               System.out.println("CVode failed with flag = " + flag);\r
+                               break;\r
+                       }\r
+                       \r
+                       if (flag == CV_SUCCESS) {\r
+                               iout++;\r
+                               tout *= TMULT;\r
+                       }\r
+                       \r
+                       if (iout == NOUT) {\r
+                               break;\r
+                       }\r
+               } // while (true)\r
+               \r
+               // Print some final statistics\r
+               PrintFinalStats(cvode_mem);\r
+               \r
+               // Check the solution error\r
+               flag = check_ans(y, reltol, abstol);\r
+               \r
+               // Free y and abstol vectors\r
+               N_VDestroy(y);\r
+               N_VDestroy(abstol);\r
+               \r
+               // Free the integrator memory\r
+               CVodeFree(cvode_mem);\r
+               \r
+               // Free the linear solver memory\r
+               SUNLinSolFree_Dense(LS);\r
+               \r
+               // Free the matrix memory\r
+               for (i = 0; i < NEQ; i++) {\r
+                       A[i] = null;\r
+               }\r
+               A = null;\r
+               return;\r
+       }\r
+       \r
        \r
 \r
        private void PrintFinalStats(CVodeMemRec cv_mem) {\r
@@ -1297,8 +1575,8 @@ public abstract class CVODES {
                System.out.println("xdotdot - 3*(1 - x^2)*xdot + x = 0, x(0) = 2, xdot(0) = 0");\r
                System.out.println("neq = " + P1_NEQ + ", reltol = " + RTOL + ", abstol = " + ATOL);\r
                \r
-               problem = cvSDirectDemo_ls_Problem_1;\r
-               f1 = cvSDirectDemo_ls_Problem_1;\r
+               problem = cvsDirectDemo_ls_Problem_1;\r
+               f1 = cvsDirectDemo_ls_Problem_1;\r
                cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);\r
                if (cvode_mem == null) {\r
                        return;\r
@@ -1473,7 +1751,7 @@ public abstract class CVODES {
                        if (flag < 0) {\r
                                return flag;\r
                        }\r
-                       Jac1 = cvSDirectDemo_ls_Problem_1;\r
+                       Jac1 = cvsDirectDemo_ls_Problem_1;\r
                        flag = CVDlsSetJacFn(cvode_mem, Jac1);\r
                                if (flag != CVDLS_SUCCESS) {\r
                                        return flag;\r
@@ -1576,28 +1854,28 @@ public abstract class CVODES {
         * @param user_data\r
         * @return\r
         */\r
-       private int fTestMode(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
+       private 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)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)  || (problem == cvsRoberts_FSA_dns)) {\r
                    ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
                    ydot[2] = 3.0E7*y[1]*y[1];\r
                    ydot[1] = -ydot[0] - ydot[2];\r
                }\r
-               else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+               else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                        ydot[0] = y[1];\r
                        ydot[1] = (ONE - y[0]*y[0]) * P1_ETA * y[1] - y[0];\r
                }\r
                return 0;\r
        }\r
        \r
-       public abstract int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data);\r
+       public abstract int f(double t, NVector yv, NVector ydotv, UserData user_data);\r
        \r
-       private int fQTestMode(double t, NVector x, NVector y, CVodeMemRec user_data) {\r
+       private int fQTestMode(double t, NVector x, NVector y, UserData user_data) {\r
                return 0;\r
        }\r
        \r
-       public abstract int fQ(double t, NVector x, NVector y, CVodeMemRec user_data);\r
+       public abstract int fQ(double t, NVector x, NVector y, UserData user_data);\r
        \r
        /**\r
         * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
@@ -1607,9 +1885,9 @@ public abstract class CVODES {
         * @param user_data\r
         * @return\r
         */\r
-       private int gTestMode(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
+       private int gTestMode(double t, NVector yv, double gout[], UserData user_data) {\r
                double y[] = yv.getData();\r
-               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_dns)) {\r
                    gout[0] = y[0] - 0.0001;\r
                    gout[1] = y[2] - 0.01;\r
                }\r
@@ -1617,7 +1895,7 @@ public abstract class CVODES {
                return 0;\r
        }\r
        \r
-       public abstract int g(double t, NVector yv, double gout[], CVodeMemRec user_data);\r
+       public abstract int g(double t, NVector yv, double gout[], UserData user_data);\r
        \r
        /**\r
         * Jacobian routine.  Compute J(t,y) = df/dy.\r
@@ -1629,9 +1907,9 @@ public abstract class CVODES {
         * @param tmp2\r
         * @return\r
         */\r
-       private int JacTestMode(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
+       private int JacTestMode(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
                double y[] = yv.getData();\r
-               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_dns)) {\r
                    J[0][0] = -0.04;\r
                    J[0][1] = 1.0E4 * y[2];\r
                    J[0][2] = 1.0E4 * y[1];\r
@@ -1644,7 +1922,7 @@ public abstract class CVODES {
                    J[2][1] = 6.0E7 * y[1];\r
                    J[2][2] = ZERO;\r
                }\r
-               else if (problem == cvSDirectDemo_ls_Problem_1) {\r
+               else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                        J[0][0] = ZERO;\r
                        J[0][1] = ONE;\r
                        J[1][0] = -TWO * P1_ETA * y[0] * y[1] - ONE;\r
@@ -1654,10 +1932,10 @@ public abstract class CVODES {
            return 0;\r
        }\r
        \r
-       public abstract int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3);\r
+       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, CVodeMemRec user_data) {\r
-               if (problem == cvsRoberts_dns_uw) {\r
+       public int ewtTestMode(NVector y, NVector w, UserData user_data) {\r
+               if ((problem == cvsRoberts_dns_uw) || (problem == cvsRoberts_FSA_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
@@ -1687,7 +1965,7 @@ public abstract class CVODES {
                return 0;\r
        }\r
        \r
-       public abstract int ewt(NVector y, NVector w, CVodeMemRec user_data);\r
+       public abstract int ewt(NVector y, NVector w, UserData user_data);\r
        \r
        \r
        \r
@@ -1708,6 +1986,11 @@ public abstract class CVODES {
                }\r
        }\r
        \r
+       public class UserData {\r
+               CVodeMemRec memRec;\r
+               double[] array;\r
+       }\r
+       \r
          \r
        public class CVodeMemRec {\r
            \r
@@ -1720,7 +2003,7 @@ public abstract class CVODES {
          //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
          int cv_f;\r
          //void *cv_user_data;         /* user pointer passed to f                      */\r
-         CVodeMemRec cv_user_data; \r
+         UserData cv_user_data = new UserData(); \r
 \r
          int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
          int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
@@ -1733,7 +2016,7 @@ public abstract class CVODES {
          //CVEwtFn cv_efun;            /* function to set ewt                           */\r
          int cv_efun_select;\r
          //void *cv_e_data;            /* user pointer passed to efun                   */\r
-         CVodeMemRec cv_e_data;\r
+         UserData cv_e_data = new UserData();\r
 \r
          /*-----------------------\r
            Quadrature Related Data \r
@@ -2162,7 +2445,7 @@ public abstract class CVODES {
          cv_mem.cv_user_efun  = false;\r
          //cv_mem.cv_efun       = null;\r
          cv_mem.cv_efun_select = -1;\r
-         cv_mem.cv_e_data     = null;\r
+         cv_mem.cv_e_data     =  new UserData();\r
          //cv_mem.cv_ehfun      = cvErrHandler;\r
          cv_mem.cv_eh_data    = cv_mem;\r
          //cv_mem.cv_errfp      = stderr;\r
@@ -2609,7 +2892,7 @@ public abstract class CVODES {
        cv_mem.cv_user_efun = false;\r
        //cv_mem.cv_efun = cvEwtSet;\r
        cv_mem.cv_efun_select = cvEwtSet_select;\r
-       cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
+       cv_mem.cv_e_data = new UserData(); /* will be set to cvode_mem in InitialSetup */\r
 \r
        return(CV_SUCCESS);\r
      }\r
@@ -2634,7 +2917,7 @@ public abstract class CVODES {
          cv_mem.cv_user_efun = true;\r
          //cv_mem.cv_efun = efun;\r
          cv_mem.cv_efun_select = efun_select;\r
-         cv_mem.cv_e_data = null; /* will be set to user_data in InitialSetup */\r
+         cv_mem.cv_e_data = new UserData(); /* will be set to user_data in InitialSetup */\r
 \r
          return(CV_SUCCESS);\r
        }\r
@@ -2903,7 +3186,7 @@ public abstract class CVODES {
       /* Initialize Jacobian-related data */\r
       cvdls_mem.jacDQ = true;\r
       cvdls_mem.jac = cvDlsDQJac;\r
-      cvdls_mem.J_data = cv_mem;\r
+      cvdls_mem.J_data.memRec = cv_mem;\r
       cvdls_mem.last_flag = CVDLS_SUCCESS;\r
 \r
       /* Initialize counters */\r
@@ -2960,7 +3243,7 @@ public abstract class CVODES {
     //CVDlsJacFn jac;       /* Jacobian routine to be called                 */\r
     int jac;\r
     //void *J_data;         /* data pointer passed to jac                    */\r
-    CVodeMemRec J_data;\r
+    UserData J_data = new UserData();\r
 \r
     double A[][];          /* A = I - gamma * df/dy                         */\r
     double savedJ[][];     /* savedJ = old Jacobian                         */\r
@@ -3006,7 +3289,7 @@ public abstract class CVODES {
     } else {\r
       cvdls_mem.jacDQ  = true;\r
       cvdls_mem.jac    = cvDlsDQJac;\r
-      cvdls_mem.J_data = cv_mem;\r
+      cvdls_mem.J_data.memRec = cv_mem;\r
     }\r
 \r
     return(CVDLS_SUCCESS);\r
@@ -3645,7 +3928,7 @@ public abstract class CVODES {
 \r
     /* Set data for efun */\r
     if (cv_mem.cv_user_efun) cv_mem.cv_e_data = cv_mem.cv_user_data;\r
-    else                      cv_mem.cv_e_data = cv_mem;\r
+    else                      cv_mem.cv_e_data.memRec = cv_mem;\r
 \r
     /* Load initial error weights */\r
     if (testMode && cv_mem.cv_itol == CV_WF) {\r
@@ -4297,11 +4580,11 @@ public abstract class CVODES {
 \r
 \r
 \r
-  private int cv_efun(NVector ycur, NVector weight, CVodeMemRec cv_mem, int cv_efun_select) {\r
+  private int cv_efun(NVector ycur, NVector weight, UserData user_data, int cv_efun_select) {\r
          int flag = 0;\r
          switch (cv_efun_select) {\r
          case cvEwtSet_select:\r
-                 flag = cvEwtSet(ycur, weight, cv_mem);  \r
+                 flag = cvEwtSet(ycur, weight, user_data.memRec);  \r
                  break;\r
          }\r
          return flag;\r
@@ -4579,7 +4862,7 @@ public abstract class CVODES {
      it has changed based on user input) */\r
   if (cvdls_mem.jacDQ) {\r
     cvdls_mem.jac    = cvDlsDQJac;\r
-    cvdls_mem.J_data = cv_mem;\r
+    cvdls_mem.J_data.memRec = cv_mem;\r
   } else {\r
     cvdls_mem.J_data = cv_mem.cv_user_data;\r
   }\r
@@ -8983,7 +9266,7 @@ else                return(snrm);
           /* Save quadrature correction in acorQ */\r
           retval = cvQuadSensRhsInternalDQ(cv_mem.cv_Ns, cv_mem.cv_tn, cv_mem.cv_y,\r
                                   cv_mem.cv_yS, cv_mem.cv_ftempQ,\r
-                                  cv_mem.cv_acorQS, cv_mem.cv_user_data,\r
+                                  cv_mem.cv_acorQS, cv_mem.cv_user_data.memRec,\r
                                   cv_mem.cv_tempv, cv_mem.cv_tempvQ);\r
           cv_mem.cv_nfQSe++;\r
           if (retval < 0) return(CV_QSRHSFUNC_FAIL);\r
@@ -10082,7 +10365,7 @@ else                return(snrm);
             cv_mem.cv_user_efun = false;\r
             //cv_mem.cv_efun = cvEwtSet;\r
             cv_mem.cv_efun_select = cvEwtSet_select;\r
-            cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
+            cv_mem.cv_e_data = new UserData(); /* will be set to cvode_mem in InitialSetup */\r
 \r
             return(CV_SUCCESS);\r
           }\r
@@ -10352,6 +10635,82 @@ else                return(snrm);
          v.own_data = false;\r
          return(v);\r
        }\r
+       \r
+       /* ----------------------------------------------------------------------------\r
+        * Function to create an array of new serial vectors. \r
+        */\r
+\r
+       private NVector[] N_VCloneVectorArray_Serial(int count, NVector w)\r
+       {\r
+         NVector vs[];\r
+         int j;\r
+\r
+         if (count <= 0) return(null);\r
+\r
+         vs = null;\r
+         vs = new NVector[count];\r
+         if(vs == null) return(null);\r
+\r
+         for (j = 0; j < count; j++) {\r
+           vs[j] = null;\r
+           vs[j] = N_VClone_Serial(w);\r
+           if (vs[j] == null) {\r
+             N_VDestroyVectorArray_Serial(vs, j-1);\r
+             return(null);\r
+           }\r
+         }\r
+\r
+         return(vs);\r
+       }\r
+       \r
+       private NVector N_VClone_Serial(NVector w)\r
+       {\r
+         NVector v;\r
+         double data[];\r
+         int length;\r
+\r
+         v = null;\r
+         v = N_VCloneEmpty_Serial(w);\r
+         if (v == null) return(null);\r
+\r
+         length = w.data.length;\r
+\r
+         /* Create data */\r
+         if (length > 0) {\r
+\r
+           /* Allocate memory */\r
+           data = null;\r
+           data = new double[length];\r
+\r
+           /* Attach data */\r
+           v.own_data = true;\r
+           v.data     = data;\r
+\r
+         }\r
+\r
+         return(v);\r
+       }\r
+       \r
+       /* ----------------------------------------------------------------------------\r
+        * Function to free an array created with N_VCloneVectorArray_Serial\r
+        */\r
+\r
+       private void N_VDestroyVectorArray_Serial(NVector vs[], int count)\r
+       {\r
+         int j;\r
+\r
+         for (j = 0; j < count; j++) {\r
+                 vs[j].data = null;\r
+                 vs[j] = null;\r
+         }\r
+\r
+         vs = null;\r
+\r
+         return;\r
+       }\r
+\r
+\r
+\r
 \r
 \r
 \r