Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 15 Mar 2018 21:13:35 +0000 (21:13 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 15 Mar 2018 21:13:35 +0000 (21:13 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15417 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 7e0392c..c5ebd0c 100644 (file)
@@ -4,6 +4,7 @@ import java.io.RandomAccessFile;
 \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.view.MipavUtil;\r
 import gov.nih.mipav.view.Preferences;\r
 \r
@@ -107,6 +108,8 @@ public abstract class CVODES {
        final int CV_NORMAL = 1;\r
        final int CV_ONE_STEP = 2;\r
        \r
+       \r
+       boolean sensi = true;\r
        //ism:   This parameter specifies the sensitivity corrector type\r
        //        to be used. In the CV_SIMULTANEOUS case, the nonlinear\r
        //        systems for states and all sensitivities are solved\r
@@ -118,7 +121,8 @@ public abstract class CVODES {
        final int CV_SIMULTANEOUS = 1;\r
        final int CV_STAGGERED = 2;\r
        final int CV_STAGGERED1 = 3;\r
-\r
+    int sensi_meth = CV_SIMULTANEOUS;\r
+    boolean err_con = true;\r
 \r
        \r
        /*\r
@@ -561,6 +565,10 @@ public abstract class CVODES {
     final int cvDlsFree_select = 1;\r
     // For cv_mem.cv_setup_select\r
     final int cvDlsSetup_select = 1;\r
+    // For cv_mem.cv_fS1_select\r
+    final int cvSensRhs1InternalDQ_select = -1;\r
+    // For cv_mem.cv_f_select\r
+    final int cvSensRhsInternalDQ_select = -1;\r
 \r
     final int cvsRoberts_dns = 1;\r
     final int cvsDirectDemo_ls_Problem_1 = 2;\r
@@ -603,6 +611,11 @@ public abstract class CVODES {
                  return 0;\r
          }\r
          \r
+         public int fS(int Ns, double t, NVector yv, NVector ydot, int is,\r
+                                NVector yS, NVector ySdot, UserData user_data, NVector tmp1, NVector tmp2) {\r
+                         return 0;  \r
+          }\r
+         \r
          public int Jac(double t, NVector yv, NVector fy, double J[][], UserData data, NVector tmp1, \r
                          NVector tmp2, NVector tmp3) {\r
                  return 0;\r
@@ -969,9 +982,9 @@ public abstract class CVODES {
                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 f = cvsRoberts_dns_uw;\r
+               int g = cvsRoberts_dns_uw;\r
+               int Jac = cvsRoberts_dns_uw;\r
                int ewt_select = cvEwtUser_select1;\r
                \r
                // Set initial conditions\r
@@ -1223,9 +1236,11 @@ public abstract class CVODES {
                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 qu;\r
+               double hu;\r
+               int f = cvsRoberts_FSA_dns;\r
+               int g = cvsRoberts_FSA_dns;\r
+               int Jac = cvsRoberts_FSA_dns;\r
                int ewt_select = cvEwtUser_select1;\r
                // User data structure\r
                double p[] = new double[]{0.04,1.0E4,3.0E7};\r
@@ -1234,7 +1249,7 @@ public abstract class CVODES {
                \r
                // Set initial conditions\r
                NVector y = new NVector();\r
-               NVector yS[];\r
+               NVector yS[] = null;\r
                NVector abstol = new NVector();\r
                double yr[] = new double[]{Y1,Y2,Y3};\r
                N_VNew_Serial(y, NEQ);\r
@@ -1254,9 +1269,9 @@ public abstract class CVODES {
                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
+               int is;\r
+               int fS = cvsRoberts_FSA_dns;            \r
                // Call CVodeCreate to create the solver memory and specify the\r
                // Backward Differentiation Formula and the use of a Newton\r
                // iteration\r
@@ -1340,78 +1355,132 @@ public abstract class CVODES {
                        if (yS == null) {\r
                                return;\r
                        }\r
+                       \r
+                       for (is = 0; is < NS; is++) {\r
+                           N_VConst_Serial(ZERO, yS[is]);      \r
+                       }\r
+                       \r
+                       // Call CVodeSensInit1 to activate forward sensitivity computations\r
+                       // and allocate internal memory for COVEDS related to sensitivity\r
+                       // calculations.  Computes the right-hand sides of the sensitivity\r
+                       // ODE, one at a time.\r
+                       flag = CVodeSensInit1(cvode_mem, NS, sensi_meth, fS, yS);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       \r
+                       // Call CVodeSensEEtolerances to estimate tolerances for sensitivity\r
+                       // variables based on the tolerances supplied for state variables and\r
+                       // the scaling factor pbar.\r
+                       flag = CVodeSensEEtolerances(cvode_mem);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       \r
+                       // Set sensitivity analysis optional inputs.\r
+                       // Call CVodeSetSensErrCon to specify the error control strategy for\r
+                       // sensitivity variables.\r
+                       cvode_mem.cv_errconS = err_con;\r
+                       \r
+                       // Call CVodeSetSensParams to specify problem parameter information for\r
+                       // sensitivity calculations.\r
+                       flag = CVodeSetSensParams(cvode_mem, null, pbar, null);\r
+                       if (flag < 0) {\r
+                               return;\r
+                       }\r
+                       \r
+                       System.out.print("Sensitivity: YES ");\r
+                       if (sensi_meth == CV_SIMULTANEOUS) {\r
+                               System.out.print("( SIMULTANEOUS + ");\r
+                       }\r
+                       else if (sensi_meth == CV_STAGGERED) {\r
+                               System.out.print("( STAGGERED + ");\r
+                       }\r
+                       else {\r
+                               System.out.print("( STAGGERED1 + ");\r
+                       }\r
+                       if (err_con) {\r
+                               System.out.println(" FULL ERROR CONTROL )");\r
+                       }\r
+                       else {\r
+                               System.out.println(" PARTIAL ERROR CONTROL )");\r
+                       }\r
                } // if (sensi)\r
+               else {\r
+                       System.out.println("Sensitivity: NO ");\r
+               }\r
                \r
-               iout = 0;\r
-               tout = T1;\r
+               // In loop over output points, call CVode, print results, test for error.\r
                \r
-               while (true) {\r
+               for (iout = 1, tout = T1; iout <= NOUT; iout++, tout *= TMULT) {\r
                        flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       if (flag < 0) {\r
+                               break;\r
+                       }\r
                        switch (iout) {\r
-                       case 0:\r
+                       case 1:\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
+                       case 2:\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
+                       case 3:\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
+                       case 4:\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
+                       case 5:\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
+                       case 6:\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
+                       case 7:\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
+                       case 8:\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
+                       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 9:\r
+                       case 910:\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
+                       case 11:\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
+                       case 12:\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("Number of integrations steps = " + cvode_mem.cv_nst);\r
+                       // Returns the order on the last successful step\r
+                   qu = cvode_mem.cv_qu;\r
+                   // Return the step size used on the last successful step\r
+                   hu = cvode_mem.cv_hu;\r
+                   System.out.println("Step order qu = " + qu);\r
+                   System.out.println("Step size hu = " + hu);\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
+                       // Call CVodeGetSens to get the sensitivity solution vector after a\r
+                       // successful return from CVode\r
+                       if (sensi) {\r
+                               flag = CVodeGetSens(cvode_mem, t, yS);\r
+                               if (flag < 0) {\r
+                                   break;      \r
+                               }\r
+                               System.out.println("Sensitivity 1 = " + yS[0].data[0] + "  " + yS[0].data[1] + "  " + yS[0].data[2]);\r
+                               System.out.println("Sensitivity 2 = " + yS[1].data[0] + "  " + yS[1].data[1] + "  " + yS[1].data[2]);\r
+                               System.out.println("Sensitivity 3 = " + yS[2].data[0] + "  " + yS[2].data[1] + "  " + yS[2].data[2]);\r
+                       } // if (sensi)\r
                        \r
-                       if (iout == NOUT) {\r
-                               break;\r
-                       }\r
-               } // while (true)\r
+               } // for (iout = 1, tout = T1; iout <= NOUT; iout++, tout *= TMULT)\r
                \r
                // Print some final statistics\r
-               PrintFinalStats(cvode_mem);\r
+               PrintFinalStats(cvode_mem, sensi);\r
                \r
                // Check the solution error\r
                flag = check_ans(y, reltol, abstol);\r
@@ -1448,6 +1517,32 @@ public abstract class CVODES {
            System.out.println("Number of calls to g (for rootfinding) = " + cv_mem.cv_nge);\r
        }\r
        \r
+       private void PrintFinalStats(CVodeMemRec cv_mem, boolean sensi) {\r
+               long nfeLS;\r
+           System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
+           System.out.println("Number of calls to f = " + cv_mem.cv_nfe);\r
+           System.out.println("Number of calls to the linear solver setup routine = " + cv_mem.cv_nsetups);\r
+           System.out.println("Number of error test failures = " + cv_mem.cv_netf[0]);\r
+           System.out.println("Number of nonlinear solver iterations = " + cv_mem.cv_nni);\r
+           System.out.println("Number of nonlinear solver convergence failures = " + cv_mem.cv_ncfn[0]);\r
+           \r
+           if (sensi) {\r
+               System.out.println("Number of calls to the sensitivity right hand side routine = " + cv_mem.cv_nfSe);   \r
+               System.out.println("Number of calls to the user f routine due to finite difference evaluations ");\r
+               System.out.println("of the sensitivity equations = " + cv_mem.cv_nfeS);\r
+            System.out.println("Number of calls made to the linear solver's setup routine due to ");\r
+            System.out.println("sensitivity computations = " + cv_mem.cv_nsetupsS);\r
+            System.out.println("Number of local error test failures for sensitivity variables = " + cv_mem.cv_netfS);\r
+            System.out.println("Total number of nonlinear iterations for sensitivity variables = " + cv_mem.cv_nniS);\r
+            System.out.println("Total number of nonlinear convergence failures for sensitivity variables = " + cv_mem.cv_ncfnS);\r
+           } // if (sensi)\r
+           \r
+           System.out.println("Number of Jacobian evaluations = " + cv_mem.cv_lmem.nje);\r
+           // the number of calls to the ODE function needed for the DQ Jacobian approximation\r
+               nfeLS = cv_mem.cv_lmem.nfeDQ;\r
+               System.out.println("Number of f evaluations in linear solver = " + nfeLS);\r
+       }\r
+       \r
        /* compare the solution at the final time 4e10s to a reference solution computed\r
           using a relative tolerance of 1e-8 and absolute tolerance of 1e-14 */\r
        private int check_ans(NVector y, double rtol, NVector atol)\r
@@ -1857,11 +1952,17 @@ public abstract class CVODES {
        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)  || (problem == cvsRoberts_FSA_dns)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\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 == cvsRoberts_FSA_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
+               ydot[1] =-ydot[0] - ydot[2];\r
+               }\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
@@ -1877,6 +1978,39 @@ public abstract class CVODES {
        \r
        public abstract int fQ(double t, NVector x, NVector y, UserData user_data);\r
        \r
+       private int fsTestMode(int Ns, double t, NVector yv, NVector ydot, int is,\r
+                       NVector yS, NVector ySdot, UserData user_data, NVector tmp1, NVector tmp2) {\r
+               double y[] = yv.data;\r
+               double p[] = user_data.array;\r
+               double s[] = yS.data;\r
+               double sd[] = ySdot.data;\r
+               if (problem == cvsRoberts_FSA_dns) {\r
+                       sd[0] = -p[0]*s[0] + p[1]*y[2]*s[1] + p[1]*y[1]*s[2];\r
+                       sd[2] = 2.0*p[2]*y[1]*s[1];\r
+                       sd[1] = -sd[0] -sd[2];\r
+                       \r
+                       switch (is) {\r
+                       case 0:\r
+                               sd[0] += -y[0];\r
+                               sd[1] += y[0];\r
+                               break;\r
+                       case 1:\r
+                               sd[0] += y[1]*y[2];\r
+                               sd[1] += -y[1]*y[2];\r
+                               break;\r
+                       case 2:\r
+                               sd[1] += -y[1]*y[1];\r
+                               sd[2] += y[1]*y[1];\r
+                               break;\r
+                       }\r
+               }\r
+               \r
+               return 0;\r
+       }\r
+       \r
+       public abstract int fS(int Ns, double t, NVector yv, NVector ydot, int is,\r
+                       NVector yS, NVector ySdot, UserData user_data, NVector tmp1, NVector tmp2);\r
+       \r
        /**\r
         * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
         * @param t\r
@@ -1909,7 +2043,7 @@ public abstract class CVODES {
         */\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) || (problem == cvsRoberts_FSA_dns)) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
                    J[0][0] = -0.04;\r
                    J[0][1] = 1.0E4 * y[2];\r
                    J[0][2] = 1.0E4 * y[1];\r
@@ -1922,6 +2056,20 @@ 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
+                       double p[] = data.array;\r
+                       J[0][0] = -p[0];\r
+                       J[0][1] = p[1]*y[2];\r
+                       J[0][2] = p[1]*y[1];\r
+                       \r
+                       J[1][0] = p[0];\r
+                       J[1][1] = -p[1]*y[2] - 2.0*p[2]*y[1];\r
+                       J[1][2] = -p[1]*y[1];\r
+                       \r
+                       J[2][0] = ZERO;\r
+                       J[2][1] = 2.0*p[2]*y[1];\r
+                       J[2][2] = ZERO;\r
+               }\r
                else if (problem == cvsDirectDemo_ls_Problem_1) {\r
                        J[0][0] = ZERO;\r
                        J[0][1] = ONE;\r
@@ -2044,9 +2192,11 @@ public abstract class CVODES {
          int cv_ism;                 /* ism = SIMULTANEOUS or STAGGERED              */\r
 \r
          //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
+         int cv_fS_select;\r
          //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
+         int cv_fS1_select;\r
          //void *cv_fS_data;           /* data pointer passed to fS                    */\r
-         CVodeMemRec cv_fS_data;\r
+         UserData cv_fS_data = new UserData();;\r
          boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
          int cv_ifS;                 /* ifS = ALLSENS or ONESENS                     */\r
 \r
@@ -2485,9 +2635,11 @@ public abstract class CVODES {
          /* Set default values for sensi. optional inputs */\r
 \r
          cv_mem.cv_sensi      = false;\r
-         //cv_mem.cv_fS_data    = mull;\r
+         cv_mem.cv_fS_data    = new UserData();\r
          //cv_mem.cv_fS         = cvSensRhsInternalDQ;\r
-         //cv_mem.cv_fS1        = cvSensRhs1InternalDQ;\r
+         cv_mem.cv_fS_select = cvSensRhsInternalDQ_select;\r
+         //cv_mem.cv_fS1         = cvSensRhs1InternalDQ;\r
+         cv_mem.cv_fS1_select = cvSensRhs1InternalDQ_select;\r
          cv_mem.cv_fSDQ       = true;\r
          cv_mem.cv_ifS        = CV_ONESENS;\r
          cv_mem.cv_DQtype     = CV_CENTERED;\r
@@ -4952,12 +5104,12 @@ public abstract class CVODES {
 \r
    if (cv_mem.cv_ifS==CV_ALLSENS) {\r
      retval = cvSensRhsInternalDQ(cv_mem.cv_Ns, time, ycur, fcur, yScur, \r
-                            fScur, cv_mem.cv_fS_data, temp1, temp2);\r
+                            fScur, cv_mem.cv_fS_data.memRec, temp1, temp2);\r
      cv_mem.cv_nfSe++;\r
    } else {\r
      for (is=0; is<cv_mem.cv_Ns; is++) {\r
        retval = cvSensRhs1InternalDQ(cv_mem.cv_Ns, time, ycur, fcur, is, yScur[is], \r
-                               fScur[is], cv_mem.cv_fS_data, temp1, temp2);\r
+                               fScur[is], cv_mem.cv_fS_data.memRec, temp1, temp2);\r
        cv_mem.cv_nfSe++;\r
        if (retval != 0) break;\r
      }\r
@@ -8882,7 +9034,7 @@ else                return(snrm);
           int retval;\r
 \r
           retval = cvSensRhs1InternalDQ(cv_mem.cv_Ns, time, ycur, fcur, is, yScur, \r
-                                  fScur, cv_mem.cv_fS_data, temp1, temp2);\r
+                                  fScur, cv_mem.cv_fS_data.memRec, temp1, temp2);\r
           cv_mem.cv_nfSe++;\r
 \r
           return(retval);\r
@@ -10710,6 +10862,484 @@ else                return(snrm);
        }\r
 \r
 \r
+       /*\r
+        * CVodeSensInit1\r
+        *\r
+        * CVodeSensInit1 allocates and initializes sensitivity related \r
+        * memory for a problem (using a sensitivity RHS function of type\r
+        * CVSensRhs1Fn). All problem specification inputs are checked for \r
+        * errors.\r
+        * The return value is CV_SUCCESS = 0 if no errors occurred, or\r
+        * a negative value otherwise.\r
+        */\r
+    // CVSensRhs1Fn fS1\r
+       private int CVodeSensInit1(CVodeMemRec cv_mem, int Ns, int ism, int fS1_select, NVector yS0[])\r
+       {\r
+         boolean allocOK;\r
+         int is;\r
+         int i;\r
+         \r
+         /* Check cvode_mem */\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Check if CVodeSensInit or CVodeSensInit1 was already called */\r
+\r
+         if (cv_mem.cv_SensMallocDone) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_SENSINIT_2);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Check if Ns is legal */\r
+\r
+         if (Ns<=0) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_BAD_NS);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         cv_mem.cv_Ns = Ns;\r
+\r
+         /* Check if ism is legal */\r
+\r
+         if ((ism!=CV_SIMULTANEOUS) && (ism!=CV_STAGGERED) && (ism!=CV_STAGGERED1)) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_BAD_ISM);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+         cv_mem.cv_ism = ism;\r
+\r
+         /* Check if yS0 is non-null */\r
+\r
+         if (yS0 == null) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_NULL_YS0);\r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         /* Store sensitivity RHS-related data */\r
+\r
+         cv_mem.cv_ifS = CV_ONESENS;\r
+         //cv_mem.cv_fS  = null;\r
+         cv_mem.cv_fS_select = -1;\r
+\r
+         if (fS1_select < 0) {\r
+\r
+           cv_mem.cv_fSDQ = true;\r
+           cv_mem.cv_fS1_select  = cvSensRhs1InternalDQ_select;\r
+           cv_mem.cv_fS_data.memRec = cv_mem;\r
+\r
+         } else {\r
+\r
+           cv_mem.cv_fSDQ = false;\r
+           cv_mem.cv_fS1_select  = fS1_select;\r
+           cv_mem.cv_fS_data = cv_mem.cv_user_data;\r
+\r
+         }\r
+\r
+         /* Allocate ncfS1, ncfnS1, and nniS1 if needed */\r
+\r
+         if (ism == CV_STAGGERED1) {\r
+           cv_mem.cv_stgr1alloc = true;\r
+           cv_mem.cv_ncfS1 = null;\r
+           cv_mem.cv_ncfS1 = new int[Ns][1];\r
+           cv_mem.cv_ncfnS1 = null;\r
+           cv_mem.cv_ncfnS1 = new long[Ns][1];\r
+           cv_mem.cv_nniS1 = null;\r
+           cv_mem.cv_nniS1 = new long[Ns];\r
+           if ( (cv_mem.cv_ncfS1 == null) ||\r
+                (cv_mem.cv_ncfnS1 == null) ||\r
+                (cv_mem.cv_nniS1 == null) ) {\r
+             cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeSensInit1",\r
+                            MSGCV_MEM_FAIL);\r
+             return(CV_MEM_FAIL);\r
+           }\r
+         } else {\r
+           cv_mem.cv_stgr1alloc = false;\r
+         }\r
+\r
+         /* Allocate the vectors (using yS0[0] as a template) */\r
+\r
+         allocOK = cvSensAllocVectors(cv_mem, yS0[0]);\r
+         if (!allocOK) {\r
+           if (cv_mem.cv_stgr1alloc) {\r
+             for (i = 0; i < cv_mem.cv_ncfS1.length; i++) {\r
+                 cv_mem.cv_ncfS1[i] = null;\r
+             }\r
+             cv_mem.cv_ncfS1 = null;\r
+             for (i = 0; i < cv_mem.cv_ncfnS1.length; i++) {\r
+                 cv_mem.cv_ncfnS1[i] = null;\r
+             }\r
+             cv_mem.cv_ncfnS1 = null;\r
+             cv_mem.cv_nniS1 = null;\r
+           }\r
+           cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeSensInit1",\r
+                          MSGCV_MEM_FAIL);\r
+           return(CV_MEM_FAIL);\r
+         }\r
+         \r
+         /*---------------------------------------------- \r
+           All error checking is complete at this point \r
+           -----------------------------------------------*/\r
+\r
+         /* Initialize znS[0] in the history array */\r
+\r
+         for (is=0; is<Ns; is++) \r
+           N_VScale_Serial(ONE, yS0[is], cv_mem.cv_znS[0][is]);\r
+\r
+         /* Initialize all sensitivity related counters */\r
+\r
+         cv_mem.cv_nfSe     = 0;\r
+         cv_mem.cv_nfeS     = 0;\r
+         cv_mem.cv_ncfnS[0]    = 0;\r
+         cv_mem.cv_netfS[0]    = 0;\r
+         cv_mem.cv_nniS     = 0;\r
+         cv_mem.cv_nsetupsS = 0;\r
+         if (ism==CV_STAGGERED1)\r
+           for (is=0; is<Ns; is++) {\r
+             cv_mem.cv_ncfnS1[is][0] = 0;\r
+             cv_mem.cv_nniS1[is] = 0;\r
+           }\r
+\r
+         /* Set default values for plist and pbar */\r
+\r
+         for (is=0; is<Ns; is++) {\r
+           cv_mem.cv_plist[is] = is;\r
+           cv_mem.cv_pbar[is] = ONE;\r
+         }\r
+\r
+         /* Sensitivities will be computed */\r
+\r
+         cv_mem.cv_sensi = true;\r
+         cv_mem.cv_SensMallocDone = true;\r
+\r
+         /* Sensitivity initialization was successfull */\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       /*\r
+        * cvSensAllocVectors\r
+        *\r
+        * Create (through duplication) N_Vectors used for sensitivity analysis, \r
+        * using the N_Vector 'tmpl' as a template.\r
+        */\r
+\r
+       private boolean cvSensAllocVectors(CVodeMemRec cv_mem, NVector tmpl) \r
+       {\r
+         int i, j;\r
+         \r
+         /* Allocate yS */\r
+         cv_mem.cv_yS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+         if (cv_mem.cv_yS == null) {\r
+           return(false);\r
+         }\r
+\r
+         /* Allocate ewtS */\r
+         cv_mem.cv_ewtS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+         if (cv_mem.cv_ewtS == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           return(false);\r
+         }\r
+         \r
+         /* Allocate acorS */\r
+         cv_mem.cv_acorS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+         if (cv_mem.cv_acorS == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+           return(false);\r
+         }\r
+         \r
+         /* Allocate tempvS */\r
+         cv_mem.cv_tempvS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+         if (cv_mem.cv_tempvS == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_acorS, cv_mem.cv_Ns);\r
+           return(false);\r
+         }\r
+           \r
+         /* Allocate ftempS */\r
+         cv_mem.cv_ftempS = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+         if (cv_mem.cv_ftempS == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_acorS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_tempvS, cv_mem.cv_Ns);\r
+           return(false);\r
+         }\r
+         \r
+         /* Allocate znS */\r
+         for (j=0; j<=cv_mem.cv_qmax; j++) {\r
+           cv_mem.cv_znS[j] = N_VCloneVectorArray_Serial(cv_mem.cv_Ns, tmpl);\r
+           if (cv_mem.cv_znS[j] == null) {\r
+             N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+             N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+             N_VDestroyVectorArray(cv_mem.cv_acorS, cv_mem.cv_Ns);\r
+             N_VDestroyVectorArray(cv_mem.cv_tempvS, cv_mem.cv_Ns);\r
+             N_VDestroyVectorArray(cv_mem.cv_ftempS, cv_mem.cv_Ns);\r
+             for (i=0; i<j; i++)\r
+               N_VDestroyVectorArray(cv_mem.cv_znS[i], cv_mem.cv_Ns);\r
+             return(false);\r
+           }\r
+         }\r
+         \r
+         /* Allocate space for pbar and plist */\r
+         cv_mem.cv_pbar = null;\r
+         cv_mem.cv_pbar = new double[cv_mem.cv_Ns];\r
+         if (cv_mem.cv_pbar == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_acorS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_tempvS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ftempS, cv_mem.cv_Ns);\r
+           for (i=0; i<=cv_mem.cv_qmax; i++)\r
+             N_VDestroyVectorArray(cv_mem.cv_znS[i], cv_mem.cv_Ns);\r
+           return(false);\r
+         }\r
+\r
+         cv_mem.cv_plist = null;\r
+         cv_mem.cv_plist = new int[cv_mem.cv_Ns];\r
+         if (cv_mem.cv_plist == null) {\r
+           N_VDestroyVectorArray(cv_mem.cv_yS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ewtS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_acorS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_tempvS, cv_mem.cv_Ns);\r
+           N_VDestroyVectorArray(cv_mem.cv_ftempS, cv_mem.cv_Ns);\r
+           for (i=0; i<=cv_mem.cv_qmax; i++)\r
+             N_VDestroyVectorArray(cv_mem.cv_znS[i], cv_mem.cv_Ns);\r
+           cv_mem.cv_pbar = null;\r
+           return(false);\r
+         }\r
+\r
+         /* Update solver workspace lengths */\r
+         cv_mem.cv_lrw += (cv_mem.cv_qmax + 6)*cv_mem.cv_Ns*cv_mem.cv_lrw1 + cv_mem.cv_Ns;\r
+         cv_mem.cv_liw += (cv_mem.cv_qmax + 6)*cv_mem.cv_Ns*cv_mem.cv_liw1 + cv_mem.cv_Ns;\r
+         \r
+         /* Store the value of qmax used here */\r
+         cv_mem.cv_qmax_allocS = cv_mem.cv_qmax;\r
+\r
+         return(true);\r
+       }\r
+\r
+       private int CVodeSensEEtolerances(CVodeMemRec cv_mem)\r
+       {\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSensEEtolerances",\r
+                          MSGCV_NO_MEM);    \r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was sensitivity initialized? */\r
+\r
+         if (cv_mem.cv_SensMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_SENS, "CVODES", "CVodeSensEEtolerances",\r
+                          MSGCV_NO_SENSI);\r
+           return(CV_NO_SENS);\r
+         } \r
+\r
+         cv_mem.cv_itolS = CV_EE;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
+       private int CVodeSetSensParams(CVodeMemRec cv_mem, double p[], double pbar[], int plist[])\r
+       {\r
+         int is, Ns;\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSetSensParams", MSGCV_NO_MEM);    \r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         /* Was sensitivity initialized? */\r
+\r
+         if (cv_mem.cv_SensMallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_SENS, "CVODES", "CVodeSetSensParams", MSGCV_NO_SENSI);\r
+           return(CV_NO_SENS);\r
+         }\r
+\r
+         Ns = cv_mem.cv_Ns;\r
+\r
+         /* Parameters */\r
+\r
+         cv_mem.cv_p = p;\r
+\r
+         /* pbar */\r
+\r
+         if (pbar != null)\r
+           for (is=0; is<Ns; is++) {\r
+             if (pbar[is] == ZERO) {\r
+               cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSetSensParams", MSGCV_BAD_PBAR);\r
+               return(CV_ILL_INPUT);\r
+             }\r
+             cv_mem.cv_pbar[is] = Math.abs(pbar[is]);\r
+           }\r
+         else\r
+           for (is=0; is<Ns; is++)\r
+             cv_mem.cv_pbar[is] = ONE;\r
+\r
+         /* plist */\r
+\r
+         if (plist != null)\r
+           for (is=0; is<Ns; is++) {\r
+             if ( plist[is] < 0 ) {\r
+               cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSetSensParams", MSGCV_BAD_PLIST);\r
+               return(CV_ILL_INPUT);\r
+             }\r
+             cv_mem.cv_plist[is] = plist[is];\r
+           }\r
+         else\r
+           for (is=0; is<Ns; is++)\r
+             cv_mem.cv_plist[is] = is;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+       \r
+       /* \r
+        * CVodeGetSens\r
+        *\r
+        * This routine extracts sensitivity solution into ySout at the\r
+        * time at which CVode returned the solution.\r
+        * This is just a wrapper that calls CVodeSensDky with k=0.\r
+        */\r
+        \r
+       private int CVodeGetSens(CVodeMemRec cv_mem, double tret[], NVector ySout[])\r
+       {\r
+         int flag;\r
+\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetSens", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         tret[0] = cv_mem.cv_tretlast;\r
+\r
+         flag = CVodeGetSensDky(cv_mem,cv_mem.cv_tretlast,0,ySout);\r
+\r
+         return(flag);\r
+       }\r
+       \r
+       /*\r
+        * CVodeGetSensDky\r
+        *\r
+        * If the user calls directly CVodeSensDky then s must be allocated\r
+        * prior to this call. When CVodeSensDky is called by \r
+        * CVodeGetSens, only ier=CV_SUCCESS, ier=CV_NO_SENS, or \r
+        * ier=CV_BAD_T are possible.\r
+        */\r
+\r
+       private int CVodeGetSensDky(CVodeMemRec cv_mem, double t, int k, NVector dkyS[])\r
+       {\r
+         int ier=CV_SUCCESS;\r
+         int is;\r
+         \r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetSensDky", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }  \r
+         \r
+         if (dkyS == null) {\r
+           cvProcessError(cv_mem, CV_BAD_DKY, "CVODES",\r
+                          "CVodeGetSensDky", MSGCV_NULL_DKYA);\r
+           return(CV_BAD_DKY);\r
+         }\r
+         \r
+         for (is=0; is<cv_mem.cv_Ns; is++) {\r
+           ier = CVodeGetSensDky1(cv_mem,t,k,is,dkyS[is]);\r
+           if (ier!=CV_SUCCESS) break;\r
+         }\r
+         \r
+         return(ier);\r
+       }\r
+           \r
+       /*\r
+        * CVodeGetSensDky1\r
+        *\r
+        * CVodeSensDky1 computes the kth derivative of the yS[is] function at\r
+        * time t, where tn-hu <= t <= tn, tn denotes the current         \r
+        * internal time reached, and hu is the last internal step size   \r
+        * successfully used by the solver. The user may request \r
+        * is=0, 1, ..., Ns-1 and k=0, 1, ..., qu, where qu is the current\r
+        * order. The derivative vector is returned in dky. This vector \r
+        * must be allocated by the caller. It is only legal to call this         \r
+        * function after a successful return from CVode with sensitivity \r
+        * computation enabled.\r
+        */\r
+\r
+       private int CVodeGetSensDky1(CVodeMemRec cv_mem, double t, int k, int is, NVector dkyS)\r
+       { \r
+         double s, c, r;\r
+         double tfuzz, tp, tn1;\r
+         int i, j;\r
+         \r
+         /* Check all inputs for legality */\r
+         \r
+         if (cv_mem == null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }  \r
+         \r
+         if(cv_mem.cv_sensi != true) {\r
+           cvProcessError(cv_mem, CV_NO_SENS, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_NO_SENSI);\r
+           return(CV_NO_SENS);\r
+         }\r
+\r
+         if (dkyS == null) {\r
+           cvProcessError(cv_mem, CV_BAD_DKY, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_NULL_DKY);\r
+           return(CV_BAD_DKY);\r
+         }\r
+         \r
+         if ((k < 0) || (k > cv_mem.cv_q)) {\r
+           cvProcessError(cv_mem, CV_BAD_K, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_BAD_K);\r
+           return(CV_BAD_K);\r
+         }\r
+         \r
+         if ((is < 0) || (is > cv_mem.cv_Ns-1)) {\r
+           cvProcessError(cv_mem, CV_BAD_IS, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_BAD_IS);\r
+           return(CV_BAD_IS);\r
+         }\r
+         \r
+         /* Allow for some slack */\r
+         tfuzz = FUZZ_FACTOR * cv_mem.cv_uround *\r
+           (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_hu));\r
+         if (cv_mem.cv_hu < ZERO) tfuzz = -tfuzz;\r
+         tp = cv_mem.cv_tn - cv_mem.cv_hu - tfuzz;\r
+         tn1 = cv_mem.cv_tn + tfuzz;\r
+         if ((t-tp)*(t-tn1) > ZERO) {\r
+           cvProcessError(cv_mem, CV_BAD_T, "CVODES", "CVodeGetSensDky1",\r
+                          MSGCV_BAD_T);\r
+           return(CV_BAD_T);\r
+         }\r
+         \r
+         /* Sum the differentiated interpolating polynomial */\r
+         \r
+         s = (t - cv_mem.cv_tn) / cv_mem.cv_h;\r
+         for (j=cv_mem.cv_q; j >= k; j--) {\r
+           c = ONE;\r
+           for (i=j; i >= j-k+1; i--) c *= i;\r
+           if (j == cv_mem.cv_q) {\r
+             N_VScale_Serial(c, cv_mem.cv_znS[cv_mem.cv_q][is], dkyS);\r
+           } else {\r
+             N_VLinearSum_Serial(c, cv_mem.cv_znS[j][is], s, dkyS, dkyS);\r
+           }\r
+         }\r
+         if (k == 0) return(CV_SUCCESS);\r
+         r = Math.pow(cv_mem.cv_h, -k);\r
+         N_VScale_Serial(r, dkyS, dkyS);\r
+         return(CV_SUCCESS);\r
+         \r
+       }\r
 \r
 \r
 \r