Started work on.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 22 Feb 2018 22:28:57 +0000 (22:28 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 22 Feb 2018 22:28:57 +0000 (22:28 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15388 ba61647d-9d00-f842-95cd-605cb4296b96

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

diff --git a/mipav/src/gov/nih/mipav/model/algorithms/CVODES.java b/mipav/src/gov/nih/mipav/model/algorithms/CVODES.java
new file mode 100644 (file)
index 0000000..dffd090
--- /dev/null
@@ -0,0 +1,1201 @@
+package gov.nih.mipav.model.algorithms;\r
+\r
+import java.io.RandomAccessFile;\r
+\r
+import gov.nih.mipav.view.MipavUtil;\r
+import gov.nih.mipav.view.Preferences;\r
+\r
+public abstract class CVODES {\r
+       // This is a port of code from CVODES, a stiff and non-stiff ODE solver, from C to Java\r
+       /*Copyright (c) 2002-2016, Lawrence Livermore National Security. \r
+       Produced at the Lawrence Livermore National Laboratory.\r
+       Written by A.C. Hindmarsh, D.R. Reynolds, R. Serban, C.S. Woodward,\r
+       S.D. Cohen, A.G. Taylor, S. Peles, L.E. Banks, and D. Shumaker.\r
+       LLNL-CODE-667205    (ARKODE)\r
+       UCRL-CODE-155951    (CVODE)\r
+       UCRL-CODE-155950    (CVODES)\r
+       UCRL-CODE-155952    (IDA)\r
+       UCRL-CODE-237203    (IDAS)\r
+       LLNL-CODE-665877    (KINSOL)\r
+       All rights reserved. \r
+\r
+       This file is part of SUNDIALS.  For details, \r
+       see http://computation.llnl.gov/projects/sundials\r
+\r
+       Redistribution and use in source and binary forms, with or without\r
+       modification, are permitted provided that the following conditions\r
+       are met:\r
+\r
+       1. Redistributions of source code must retain the above copyright\r
+       notice, this list of conditions and the disclaimer below.\r
+\r
+       2. Redistributions in binary form must reproduce the above copyright\r
+       notice, this list of conditions and the disclaimer (as noted below)\r
+       in the documentation and/or other materials provided with the\r
+       distribution.\r
+\r
+       3. Neither the name of the LLNS/LLNL nor the names of its contributors\r
+       may be used to endorse or promote products derived from this software\r
+       without specific prior written permission.\r
+\r
+       THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \r
+       "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT \r
+       LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS \r
+       FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL \r
+       LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF \r
+       ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, \r
+       SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED \r
+       TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, \r
+       DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY \r
+       THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
+       (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
+       OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
+\r
+       Additional BSD Notice\r
+       ---------------------\r
+       1. This notice is required to be provided under our contract with\r
+       the U.S. Department of Energy (DOE). This work was produced at\r
+       Lawrence Livermore National Laboratory under Contract \r
+       No. DE-AC52-07NA27344 with the DOE.\r
+\r
+       2. Neither the United States Government nor Lawrence Livermore \r
+       National Security, LLC nor any of their employees, makes any warranty, \r
+       express or implied, or assumes any liability or responsibility for the\r
+       accuracy, completeness, or usefulness of any */\r
+       \r
+       /* Basic CVODES constants */\r
+\r
+       final int ADAMS_Q_MAX = 12;      /* max value of q for lmm == ADAMS    */\r
+       final int BDF_Q_MAX = 5;      /* max value of q for lmm == BDF      */\r
+       final int  Q_MAX = ADAMS_Q_MAX;  /* max value of q for either lmm      */\r
+       final int L_MAX  = (Q_MAX+1);    /* max value of L for either lmm      */\r
+       final int NUM_TESTS = 5;      /* number of error test quantities    */\r
+       \r
+    double DBL_EPSILON;\r
+    double UNIT_ROUNDOFF;\r
+\r
+\r
+\r
+        // lmm:   The user of the CVODES package specifies whether to use\r
+        // the CV_ADAMS or CV_BDF (backward differentiation formula)\r
+        // linear multistep method. The BDF method is recommended\r
+        // for stiff problems, and the CV_ADAMS method is recommended\r
+        // for nonstiff problems.\r
+       final int CV_ADAMS = 1;\r
+       final int CV_BDF = 2;\r
+       \r
+       //iter:  At each internal time step, a nonlinear equation must\r
+       //        be solved. The user can specify either CV_FUNCTIONAL\r
+       //        iteration, which does not require linear algebra, or a\r
+       //        CV_NEWTON iteration, which requires the solution of linear\r
+       //        systems. In the CV_NEWTON case, the user also specifies a\r
+       //        CVODE linear solver. CV_NEWTON is recommended in case of\r
+       //        stiff problems.\r
+       final int CV_FUNCTIONAL = 1;\r
+       final int CV_NEWTON = 2;\r
+       \r
+       /*\r
+        * Control constants for type of sensitivity RHS\r
+        * ---------------------------------------------\r
+        */\r
+\r
+       final int CV_ONESENS = 1;\r
+       final int CV_ALLSENS = 2;\r
+\r
+       \r
+       /*\r
+        * Control constants for tolerances\r
+        * --------------------------------\r
+        */\r
+\r
+       final int CV_NN = 0;\r
+       final int CV_SS = 1;\r
+       final int CV_SV = 2;\r
+       final int CV_WF = 3;\r
+       final int CV_EE = 4;\r
+       \r
+       /* DQtype */\r
+       final int CV_CENTERED = 1;\r
+       final int CV_FORWARD = 2;\r
+       \r
+       /* \r
+        * ----------------------------------------\r
+        * CVODES return flags\r
+        * ----------------------------------------\r
+        */\r
+\r
+       final int  CV_SUCCESS = 0;\r
+       /*#define CV_TSTOP_RETURN          1\r
+       #define CV_ROOT_RETURN           2\r
+\r
+       #define CV_WARNING              99\r
+\r
+       #define CV_TOO_MUCH_WORK        -1\r
+       #define CV_TOO_MUCH_ACC         -2\r
+       #define CV_ERR_FAILURE          -3\r
+       #define CV_CONV_FAILURE         -4\r
+\r
+       #define CV_LINIT_FAIL           -5\r
+       #define CV_LSETUP_FAIL          -6\r
+       #define CV_LSOLVE_FAIL          -7\r
+       #define CV_RHSFUNC_FAIL         -8\r
+       #define CV_FIRST_RHSFUNC_ERR    -9\r
+       #define CV_REPTD_RHSFUNC_ERR    -10\r
+       #define CV_UNREC_RHSFUNC_ERR    -11\r
+       #define CV_RTFUNC_FAIL          -12\r
+\r
+       #define CV_MEM_FAIL             -20\r
+       #define CV_MEM_NULL             -21\r
+       #define CV_ILL_INPUT            -22\r
+       #define CV_NO_MALLOC            -23\r
+       #define CV_BAD_K                -24\r
+       #define CV_BAD_T                -25\r
+       #define CV_BAD_DKY              -26\r
+       #define CV_TOO_CLOSE            -27\r
+\r
+       #define CV_NO_QUAD              -30\r
+       #define CV_QRHSFUNC_FAIL        -31\r
+       #define CV_FIRST_QRHSFUNC_ERR   -32\r
+       #define CV_REPTD_QRHSFUNC_ERR   -33\r
+       #define CV_UNREC_QRHSFUNC_ERR   -34\r
+\r
+       #define CV_NO_SENS              -40\r
+       #define CV_SRHSFUNC_FAIL        -41\r
+       #define CV_FIRST_SRHSFUNC_ERR   -42\r
+       #define CV_REPTD_SRHSFUNC_ERR   -43\r
+       #define CV_UNREC_SRHSFUNC_ERR   -44\r
+\r
+       #define CV_BAD_IS               -45\r
+\r
+       #define CV_NO_QUADSENS          -50\r
+       #define CV_QSRHSFUNC_FAIL       -51\r
+       #define CV_FIRST_QSRHSFUNC_ERR  -52\r
+       #define CV_REPTD_QSRHSFUNC_ERR  -53\r
+       #define CV_UNREC_QSRHSFUNC_ERR  -54*/\r
+\r
+       final int CV_MEM_NULL = -21;\r
+       final int CV_ILL_INPUT = -22;\r
+\r
+\r
+\r
+       final double HMIN_DEFAULT = 0.0;    /* hmin default value     */\r
+       final double HMAX_INV_DEFAULT = 0.0;    /* hmax_inv default value */\r
+       final int MXHNIL_DEFAULT =  10;             /* mxhnil default value   */\r
+       final long MXSTEP_DEFAULT = 500L;            /* mxstep default value   */\r
+       final int NLS_MAXCOR = 3;\r
+       final int MXNCF = 10;\r
+       final int MXNEF = 7;\r
+       final double CORTES = 0.1;\r
+       final double ZERO = 0.0;\r
+       \r
+       final String MSG_TIME = "t = %lg";\r
+       final String MSG_TIME_H = "t = %lg and h = %lg";\r
+       final String MSG_TIME_INT = "t = %lg is not between tcur - hu = %lg and tcur = %lg.";\r
+       final String MSG_TIME_TOUT = "tout = %lg";\r
+       final String MSG_TIME_TSTOP = "tstop = %lg";\r
+\r
+       \r
+       /* Initialization and I/O error messages */\r
+\r
+       final String MSGCV_NO_MEM = "cvode_mem = NULL illegal.";\r
+       final String MSGCV_CVMEM_FAIL = "Allocation of cvode_mem failed.";\r
+       final String MSGCV_MEM_FAIL = "A memory request failed.";\r
+       final String MSGCV_BAD_LMM = "Illegal value for lmm. The legal values are CV_ADAMS and CV_BDF.";\r
+       final String MSGCV_BAD_ITER = "Illegal value for iter. The legal values are CV_FUNCTIONAL and CV_NEWTON.";\r
+       final String MSGCV_NO_MALLOC = "Attempt to call before CVodeInit.";\r
+       final String MSGCV_NEG_MAXORD = "maxord <= 0 illegal.";\r
+       final String MSGCV_BAD_MAXORD = "Illegal attempt to increase maximum method order.";\r
+       final String MSGCV_SET_SLDET = "Attempt to use stability limit detection with the CV_ADAMS method illegal.";\r
+       final String MSGCV_NEG_HMIN = "hmin < 0 illegal.";\r
+       final String MSGCV_NEG_HMAX = "hmax < 0 illegal.";\r
+       final String MSGCV_BAD_HMIN_HMAX = "Inconsistent step size limits: hmin > hmax.";\r
+       final String MSGCV_BAD_RELTOL = "reltol < 0 illegal.";\r
+       final String MSGCV_BAD_ABSTOL = "abstol has negative component(s) (illegal).";\r
+       final String MSGCV_NULL_ABSTOL = "abstol = NULL illegal.";\r
+       final String MSGCV_NULL_Y0 = "y0 = NULL illegal.";\r
+       final String MSGCV_NULL_F = "f = NULL illegal.";\r
+       final String MSGCV_NULL_G = "g = NULL illegal.";\r
+       final String MSGCV_BAD_NVECTOR = "A required vector operation is not implemented.";\r
+       final String MSGCV_BAD_K = "Illegal value for k.";\r
+       final String MSGCV_NULL_DKY = "dky = NULL illegal.";\r
+       final String MSGCV_BAD_T = "Illegal value for t.";\r
+       final String MSGCV_NO_ROOT = "Rootfinding was not initialized.";\r
+\r
+       final String MSGCV_NO_QUAD  = "Quadrature integration not activated.";\r
+       final String MSGCV_BAD_ITOLQ = "Illegal value for itolQ. The legal values are CV_SS and CV_SV.";\r
+       final String MSGCV_NULL_ABSTOLQ = "abstolQ = NULL illegal.";\r
+       final String MSGCV_BAD_RELTOLQ = "reltolQ < 0 illegal.";\r
+       final String MSGCV_BAD_ABSTOLQ = "abstolQ has negative component(s) (illegal).";  \r
+\r
+       final String MSGCV_SENSINIT_2 = "Sensitivity analysis already initialized.";\r
+       final String MSGCV_NO_SENSI  = "Forward sensitivity analysis not activated.";\r
+       final String MSGCV_BAD_ITOLS = "Illegal value for itolS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
+       final String MSGCV_NULL_ABSTOLS = "abstolS = NULL illegal.";\r
+       final String MSGCV_BAD_RELTOLS = "reltolS < 0 illegal.";\r
+       final String MSGCV_BAD_ABSTOLS = "abstolS has negative component(s) (illegal).";  \r
+       final String MSGCV_BAD_PBAR = "pbar has zero component(s) (illegal).";\r
+       final String MSGCV_BAD_PLIST = "plist has negative component(s) (illegal).";\r
+       final String MSGCV_BAD_NS = "NS <= 0 illegal.";\r
+       final String MSGCV_NULL_YS0 = "yS0 = NULL illegal.";\r
+       final String MSGCV_BAD_ISM = "Illegal value for ism. Legal values are: CV_SIMULTANEOUS, CV_STAGGERED and CV_STAGGERED1.";\r
+       final String MSGCV_BAD_IFS = "Illegal value for ifS. Legal values are: CV_ALLSENS and CV_ONESENS.";\r
+       final String MSGCV_BAD_ISM_IFS = "Illegal ism = CV_STAGGERED1 for CVodeSensInit.";\r
+       final String MSGCV_BAD_IS = "Illegal value for is.";\r
+       final String MSGCV_NULL_DKYA = "dkyA = NULL illegal.";\r
+       final String MSGCV_BAD_DQTYPE = "Illegal value for DQtype. Legal values are: CV_CENTERED and CV_FORWARD.";\r
+       final String MSGCV_BAD_DQRHO = "DQrhomax < 0 illegal.";\r
+\r
+       final String MSGCV_BAD_ITOLQS = "Illegal value for itolQS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
+       final String MSGCV_NULL_ABSTOLQS = "abstolQS = NULL illegal.";\r
+       final String MSGCV_BAD_RELTOLQS = "reltolQS < 0 illegal.";\r
+       final String MSGCV_BAD_ABSTOLQS = "abstolQS has negative component(s) (illegal).";  \r
+       final String MSGCV_NO_QUADSENSI = "Forward sensitivity analysis for quadrature variables not activated.";\r
+       final String MSGCV_NULL_YQS0 = "yQS0 = NULL illegal.";\r
+\r
+       /* CVode Error Messages */\r
+\r
+       final String MSGCV_NO_TOL = "No integration tolerances have been specified.";\r
+       final String MSGCV_LSOLVE_NULL = "The linear solver's solve routine is NULL.";\r
+       final String MSGCV_YOUT_NULL = "yout = NULL illegal.";\r
+       final String MSGCV_TRET_NULL = "tret = NULL illegal.";\r
+       final String MSGCV_BAD_EWT = "Initial ewt has component(s) equal to zero (illegal).";\r
+       final String MSGCV_EWT_NOW_BAD = "At " + MSG_TIME + ", a component of ewt has become <= 0.";\r
+       final String MSGCV_BAD_ITASK = "Illegal value for itask.";\r
+       final String MSGCV_BAD_H0 = "h0 and tout - t0 inconsistent.";\r
+       final String MSGCV_BAD_TOUT = "Trouble interpolating at " + MSG_TIME_TOUT + ". tout too far back in direction of integration";\r
+       final String MSGCV_EWT_FAIL = "The user-provide EwtSet function failed.";\r
+       final String MSGCV_EWT_NOW_FAIL = "At " + MSG_TIME + ", the user-provide EwtSet function failed.";\r
+       final String MSGCV_LINIT_FAIL = "The linear solver's init routine failed.";\r
+       final String MSGCV_HNIL_DONE = "The above warning has been issued mxhnil times and will not be issued again for this problem.";\r
+       final String MSGCV_TOO_CLOSE = "tout too close to t0 to start integration.";\r
+       final String MSGCV_MAX_STEPS = "At " + MSG_TIME + ", mxstep steps taken before reaching tout.";\r
+       final String MSGCV_TOO_MUCH_ACC = "At " + MSG_TIME + ", too much accuracy requested.";\r
+       final String MSGCV_HNIL = "Internal " + MSG_TIME_H + " are such that t + h = t on the next step. The solver will continue anyway.";\r
+       final String MSGCV_ERR_FAILS = "At " + MSG_TIME_H + ", the error test failed repeatedly or with |h| = hmin.";\r
+       final String MSGCV_CONV_FAILS = "At " + MSG_TIME_H + ", the corrector convergence test failed repeatedly or with |h| = hmin.";\r
+       final String MSGCV_SETUP_FAILED = "At " + MSG_TIME + ", the setup routine failed in an unrecoverable manner.";\r
+       final String MSGCV_SOLVE_FAILED = "At " + MSG_TIME + ", the solve routine failed in an unrecoverable manner.";\r
+       final String MSGCV_RHSFUNC_FAILED = "At " + MSG_TIME + ", the right-hand side routine failed in an unrecoverable manner.";\r
+       final String MSGCV_RHSFUNC_UNREC = "At " + MSG_TIME + ", the right-hand side failed in a recoverable manner, but no recovery is possible.";\r
+       final String MSGCV_RHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable right-hand side function errors.";\r
+       final String MSGCV_RHSFUNC_FIRST = "The right-hand side routine failed at the first call.";\r
+       final String MSGCV_RTFUNC_FAILED = "At " + MSG_TIME + ", the rootfinding routine failed in an unrecoverable manner.";\r
+       final String MSGCV_CLOSE_ROOTS = "Root found at and very near " + MSG_TIME + ".";\r
+       final String MSGCV_BAD_TSTOP = "The value " + MSG_TIME_TSTOP + " is behind current " + MSG_TIME + " in the direction of integration.";\r
+       final String MSGCV_INACTIVE_ROOTS = "At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.";\r
+\r
+       final String MSGCV_NO_TOLQ = "No integration tolerances for quadrature variables have been specified.";\r
+       final String MSGCV_BAD_EWTQ = "Initial ewtQ has component(s) equal to zero (illegal).";\r
+       final String MSGCV_EWTQ_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQ has become <= 0.";\r
+       final String MSGCV_QRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature right-hand side routine failed in an unrecoverable manner.";\r
+       final String MSGCV_QRHSFUNC_UNREC = "At " + MSG_TIME + ", the quadrature right-hand side failed in a recoverable manner, but no recovery is possible.";\r
+       final String MSGCV_QRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature right-hand side function errors.";\r
+       final String MSGCV_QRHSFUNC_FIRST = "The quadrature right-hand side routine failed at the first call.";\r
+\r
+       final String MSGCV_NO_TOLS = "No integration tolerances for sensitivity variables have been specified.";\r
+       final String MSGCV_NULL_P = "p = NULL when using internal DQ for sensitivity RHS illegal.";\r
+       final String MSGCV_BAD_EWTS = "Initial ewtS has component(s) equal to zero (illegal).";\r
+       final String MSGCV_EWTS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtS has become <= 0.";\r
+       final String MSGCV_SRHSFUNC_FAILED = "At " + MSG_TIME + ", the sensitivity right-hand side routine failed in an unrecoverable manner.";\r
+       final String MSGCV_SRHSFUNC_UNREC = "At " + MSG_TIME + ", the sensitivity right-hand side failed in a recoverable manner, but no recovery is possible.";\r
+       final String MSGCV_SRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable sensitivity right-hand side function errors.";\r
+       final String MSGCV_SRHSFUNC_FIRST = "The sensitivity right-hand side routine failed at the first call.";\r
+\r
+       final String MSGCV_NULL_FQ = "CVODES is expected to use DQ to evaluate the RHS of quad. sensi., but quadratures were not initialized.";\r
+       final String MSGCV_NO_TOLQS = "No integration tolerances for quadrature sensitivity variables have been specified.";\r
+       final String MSGCV_BAD_EWTQS = "Initial ewtQS has component(s) equal to zero (illegal).";\r
+       final String MSGCV_EWTQS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQS has become <= 0.";\r
+       final String MSGCV_QSRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature sensitivity right-hand side routine failed in an unrecoverable manner.";\r
+       final String MSGCV_QSRHSFUNC_UNREC = "At " + MSG_TIME + ", the quadrature sensitivity right-hand side failed in a recoverable manner, but no recovery is possible.";\r
+       final String MSGCV_QSRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature sensitivity right-hand side function errors.";\r
+       final String MSGCV_QSRHSFUNC_FIRST = "The quadrature sensitivity right-hand side routine failed at the first call.";\r
+\r
+       /* \r
+        * =================================================================\r
+        *   C V O D E A    E R R O R    M E S S A G E S\r
+        * =================================================================\r
+        */\r
+\r
+       final String MSGCV_NO_ADJ = "Illegal attempt to call before calling CVodeAdjMalloc.";\r
+       final String MSGCV_BAD_STEPS = "Steps nonpositive illegal.";\r
+       final String MSGCV_BAD_INTERP = "Illegal value for interp.";\r
+       final String MSGCV_BAD_WHICH  = "Illegal value for which.";\r
+       final String MSGCV_NO_BCK = "No backward problems have been defined yet.";\r
+       final String MSGCV_NO_FWD = "Illegal attempt to call before calling CVodeF.";\r
+       final String MSGCV_BAD_TB0 = "The initial time tB0 for problem %d is outside the interval over which the forward problem was solved.";\r
+       final String MSGCV_BAD_SENSI = "At least one backward problem requires sensitivities, but they were not stored for interpolation.";\r
+       final String MSGCV_BAD_ITASKB = "Illegal value for itaskB. Legal values are CV_NORMAL and CV_ONE_STEP.";\r
+       final String MSGCV_BAD_TBOUT  = "The final time tBout is outside the interval over which the forward problem was solved.";\r
+       final String MSGCV_BACK_ERROR  = "Error occured while integrating backward problem # %d"; \r
+       final String MSGCV_BAD_TINTERP = "Bad t = %g for interpolation.";\r
+       final String MSGCV_WRONG_INTERP = "This function cannot be called for the specified interp type.";\r
+\r
+\r
+    final int cvsRoberts_dns = 1;\r
+    int problem = cvsRoberts_dns;\r
+       \r
+       \r
+       public CVODES() {\r
+               // eps returns the distance from 1.0 to the next larger double-precision\r
+               // number, that is, eps = 2^-52.\r
+               // epsilon = D1MACH(4)\r
+               // Machine epsilon is the smallest positive epsilon such that\r
+               // (1.0 + epsilon) != 1.0.\r
+               // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
+               // epsilon = 2.2204460e-16\r
+               // epsilon is called the largest relative spacing\r
+               DBL_EPSILON = 1.0;\r
+               double neweps = 1.0;\r
+\r
+               while (true) {\r
+\r
+                       if (1.0 == (1.0 + neweps)) {\r
+                               break;\r
+                       } else {\r
+                               DBL_EPSILON = neweps;\r
+                               neweps = neweps / 2.0;\r
+                       }\r
+               } // while(true)\r
+               \r
+           UNIT_ROUNDOFF = DBL_EPSILON;\r
+           if (problem == cvsRoberts_dns) {\r
+               runcvsRoberts_dns();\r
+           }\r
+           return;\r
+       }\r
+       \r
+       private void runcvsRoberts_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 CVODE. The problem is from\r
+                * chemical kinetics, and consists of the following three rate\r
+                * equations:         \r
+                *    dy1/dt = -.04*y1 + 1.e4*y2*y3\r
+                *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2\r
+                *    dy3/dt = 3.e7*(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 problem is stiff.\r
+                * While integrating the system, we also use the rootfinding\r
+                * feature to find the points at which y1 = 1e-4 or at which\r
+                * y3 = 0.01. This program solves the problem with the BDF method,\r
+                * Newton iteration with the SUNDENSE dense linear solver, and a\r
+                * user-supplied Jacobian routine.\r
+                * It uses a scalar relative tolerance and a vector absolute\r
+                * tolerance. Output is printed in decades from t = .4 to t = 4.e10.\r
+                * Run statistics (optional outputs) are printed at the end.\r
+                * -----------------------------------------------------------------*/\r
+               \r
+               /** Problem Constants */\r
+               final int NEQ = 3; // Number of equations\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 ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+               final double ATOL2 = 1.0E-14;\r
+               final double ATOL3 = 1.0E-6;\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
+               \r
+               // Set initial conditions\r
+               double y[] = new double[]{Y1,Y2,Y3};\r
+               double reltol = RTOL; // Set the scalar relative tolerance\r
+               // Set the vector absolute tolerance\r
+               double abstol[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               CVodeMemRec cvode_mem;\r
+               int flag;\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
+               flag = CVodeInit(cvode_mem, T0, y);\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
+               \r
+       }\r
+       \r
+       private int f(double t, double y[], double ydot[]) {\r
+               if (problem == 1) {\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
+               return 0;\r
+       }\r
+       \r
+       // Types: struct CVodeMemRec, CVodeMem\r
+       // -----------------------------------------------------------------\r
+       // The type CVodeMem is type pointer to struct CVodeMemRec.\r
+       // This structure contains fields to keep track of problem state.\r
+       \r
+         \r
+       private class CVodeMemRec {\r
+           \r
+         double cv_uround;         /* machine unit roundoff                         */   \r
+\r
+         /*-------------------------- \r
+           Problem Specification Data \r
+           --------------------------*/\r
+\r
+         //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
+         //void *cv_user_data;         /* user pointer passed to f                      */\r
+\r
+         int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
+         int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
+\r
+         int cv_itol;                /* itol = CV_SS, CV_SV, or CV_WF, or CV_NN       */\r
+         double cv_reltol;         /* relative tolerance                            */\r
+         double cv_Sabstol;        /* scalar absolute tolerance                     */\r
+         double cv_Vabstol[];        /* vector absolute tolerance                     */\r
+         boolean cv_user_efun;   /* SUNTRUE if user sets efun                     */\r
+         //CVEwtFn cv_efun;            /* function to set ewt                           */\r
+         //void *cv_e_data;            /* user pointer passed to efun                   */\r
+\r
+         /*-----------------------\r
+           Quadrature Related Data \r
+           -----------------------*/\r
+\r
+         boolean cv_quadr;       /* SUNTRUE if integrating quadratures            */\r
+\r
+         //CVQuadRhsFn cv_fQ;          /* q' = fQ(t, y(t))                              */\r
+\r
+         boolean cv_errconQ;     /* SUNTRUE if quadrs. are included in error test */\r
+\r
+         int cv_itolQ;               /* itolQ = CV_SS or CV_SV                        */\r
+         double cv_reltolQ;        /* relative tolerance for quadratures            */\r
+         double cv_SabstolQ;       /* scalar absolute tolerance for quadratures     */\r
+         double cv_VabstolQ[];       /* vector absolute tolerance for quadratures     */\r
+\r
+         /*------------------------\r
+           Sensitivity Related Data \r
+           ------------------------*/\r
+\r
+         boolean cv_sensi;       /* SUNTRUE if computing sensitivities           */\r
+\r
+         int cv_Ns;                  /* Number of sensitivities                      */\r
+\r
+         int cv_ism;                 /* ism = SIMULTANEOUS or STAGGERED              */\r
+\r
+         //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
+         //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
+         //void *cv_fS_data;           /* data pointer passed to fS                    */\r
+         boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
+         int cv_ifS;                 /* ifS = ALLSENS or ONESENS                     */\r
+\r
+         double cv_p[];             /* parameters in f(t,y,p)                       */\r
+         double cv_pbar[];          /* scale factors for parameters                 */\r
+         int cv_plist[];              /* list of sensitivities                        */\r
+         int cv_DQtype;              /* central/forward finite differences           */\r
+         double cv_DQrhomax;       /* cut-off value for separate/simultaneous FD   */\r
+\r
+         boolean cv_errconS;     /* SUNTRUE if yS are considered in err. control */\r
+\r
+         int cv_itolS;\r
+         double cv_reltolS;        /* relative tolerance for sensitivities         */\r
+         double cv_SabstolS[];      /* scalar absolute tolerances for sensi.        */\r
+         double cv_VabstolS[];      /* vector absolute tolerances for sensi.        */\r
+\r
+         /*-----------------------------------\r
+           Quadrature Sensitivity Related Data \r
+           -----------------------------------*/\r
+\r
+         boolean cv_quadr_sensi; /* SUNTRUE if computing sensitivties of quadrs. */\r
+\r
+         //CVQuadSensRhsFn cv_fQS;     /* fQS = (dfQ/dy)*yS + (dfQ/dp)                 */\r
+         //void *cv_fQS_data;          /* data pointer passed to fQS                   */\r
+         boolean cv_fQSDQ;       /* SUNTRUE if using internal DQ functions       */\r
+\r
+         boolean cv_errconQS;    /* SUNTRUE if yQS are considered in err. con.   */\r
+\r
+         int cv_itolQS;\r
+         double cv_reltolQS;       /* relative tolerance for yQS                   */\r
+         double cv_SabstolQS[];     /* scalar absolute tolerances for yQS           */\r
+         double cv_VabstolQS[];     /* vector absolute tolerances for yQS           */\r
+\r
+         /*-----------------------\r
+           Nordsieck History Array \r
+           -----------------------*/\r
+\r
+         double cv_zn[] = new double[L_MAX];      /* Nordsieck array, of size N x (q+1).\r
+                                        zn[j] is a vector of length N (j=0,...,q)\r
+                                        zn[j] = [1/factorial(j)] * h^j * \r
+                                        (jth derivative of the interpolating poly.)  */\r
+\r
+         /*-------------------\r
+           Vectors of length N \r
+           -------------------*/\r
+\r
+         double cv_ewt[];            /* error weight vector                          */\r
+         double cv_y[];              /* y is used as temporary storage by the solver.\r
+                                        The memory is provided by the user to CVode \r
+                                        where the vector is named yout.              */\r
+         double cv_acor[];           /* In the context of the solution of the\r
+                                        nonlinear equation, acor = y_n(m) - y_n(0).\r
+                                        On return, this vector is scaled to give\r
+                                        the estimated local error in y.              */\r
+         double cv_tempv[];          /* temporary storage vector                     */\r
+         double cv_ftemp[];          /* temporary storage vector                     */\r
+         \r
+         /*--------------------------\r
+           Quadrature Related Vectors \r
+           --------------------------*/\r
+\r
+         double cv_znQ[] = new double[L_MAX];     /* Nordsieck arrays for quadratures             */\r
+         double cv_ewtQ[];           /* error weight vector for quadratures          */\r
+         double cv_yQ[];             /* Unlike y, yQ is not allocated by the user    */\r
+         double cv_acorQ[];          /* acorQ = yQ_n(m) - yQ_n(0)                    */\r
+         double cv_tempvQ[];         /* temporary storage vector (~ tempv)           */\r
+\r
+         /*---------------------------\r
+           Sensitivity Related Vectors \r
+           ---------------------------*/\r
+\r
+         double cv_znS[][] = new double[L_MAX][];    /* Nordsieck arrays for sensitivities           */\r
+         double cv_ewtS[][];          /* error weight vectors for sensitivities       */\r
+         double cv_yS[][];            /* yS=yS0 (allocated by the user)               */\r
+         double cv_acorS[][];         /* acorS = yS_n(m) - yS_n(0)                    */\r
+         double cv_tempvS[][];        /* temporary storage vector (~ tempv)           */\r
+         double cv_ftempS[][];        /* temporary storage vector (~ ftemp)           */\r
+\r
+         boolean cv_stgr1alloc;  /* Did we allocate ncfS1, ncfnS1, and nniS1?    */\r
+\r
+         /*--------------------------------------\r
+           Quadrature Sensitivity Related Vectors \r
+           --------------------------------------*/\r
+\r
+         double cv_znQS[][] = new double[L_MAX][];   /* Nordsieck arrays for quadr. sensitivities    */\r
+         double cv_ewtQS[][];         /* error weight vectors for sensitivities       */\r
+         double cv_yQS[][];           /* Unlike yS, yQS is not allocated by the user  */\r
+         double cv_acorQS[][];        /* acorQS = yQS_n(m) - yQS_n(0)                 */\r
+         double cv_tempvQS[][];       /* temporary storage vector (~ tempv)           */\r
+         double cv_ftempQ[];         /* temporary storage vector (~ ftemp)           */\r
+         \r
+         /*-----------------\r
+           Tstop information\r
+           -----------------*/\r
+\r
+         boolean cv_tstopset;\r
+         double cv_tstop;\r
+\r
+         /*---------\r
+           Step Data \r
+           ---------*/\r
+\r
+         int cv_q;                    /* current order                               */\r
+         int cv_qprime;               /* order to be used on the next step\r
+                                       * qprime = q-1, q, or q+1                     */\r
+         int cv_next_q;               /* order to be used on the next step           */\r
+         int cv_qwait;                /* number of internal steps to wait before\r
+                                       * considering a change in q                   */\r
+         int cv_L;                    /* L = q + 1                                   */\r
+\r
+         double cv_hin;\r
+         double cv_h;               /* current step size                           */\r
+         double cv_hprime;          /* step size to be used on the next step       */ \r
+         double cv_next_h;          /* step size to be used on the next step       */ \r
+         double cv_eta;             /* eta = hprime / h                            */\r
+         double cv_hscale;          /* value of h used in zn                       */\r
+         double cv_tn;              /* current internal value of t                 */\r
+         double cv_tretlast;        /* last value of t returned                    */\r
+\r
+         double cv_tau[] = new double[L_MAX+1];    /* array of previous q+1 successful step\r
+                                       * sizes indexed from 1 to q+1                 */\r
+         double cv_tq[] =new double[NUM_TESTS+1]; /* array of test quantities indexed from\r
+                                       * 1 to NUM_TESTS(=5)                          */\r
+         double cv_l[] = new double[L_MAX];        /* coefficients of l(x) (degree q poly)        */\r
+\r
+         double cv_rl1;             /* the scalar 1/l[1]                           */\r
+         double cv_gamma;           /* gamma = h * rl1                             */\r
+         double cv_gammap;          /* gamma at the last setup call                */\r
+         double cv_gamrat;          /* gamma / gammap                              */\r
+\r
+         double cv_crate;           /* est. corrector conv. rate in Nls            */\r
+         double cv_crateS;          /* est. corrector conv. rate in NlsStgr        */\r
+         double cv_acnrm;           /* | acor |                                    */\r
+         double cv_acnrmQ;          /* | acorQ |                                   */\r
+         double cv_acnrmS;          /* | acorS |                                   */\r
+         double cv_acnrmQS;         /* | acorQS |                                  */\r
+         double cv_nlscoef;         /* coeficient in nonlinear convergence test    */\r
+         int  cv_mnewt;               /* Newton iteration counter                    */\r
+         int  cv_ncfS1[];              /* Array of Ns local counters for conv.  \r
+                                       * failures (used in CVStep for STAGGERED1)    */\r
+\r
+         /*------\r
+           Limits \r
+           ------*/\r
+\r
+         int cv_qmax;             /* q <= qmax                                       */\r
+         long cv_mxstep;      /* maximum number of internal steps for one \r
+                                     user call                                       */\r
+         int cv_maxcor;           /* maximum number of corrector iterations for \r
+                                     the solution of the nonlinear equation          */\r
+         int cv_maxcorS;\r
+         int cv_mxhnil;           /* max. number of warning messages issued to the\r
+                                     user that t + h == t for the next internal step */\r
+         int cv_maxnef;           /* maximum number of error test failures           */\r
+         int cv_maxncf;           /* maximum number of nonlinear conv. failures      */\r
+         \r
+         double cv_hmin;        /* |h| >= hmin                                     */\r
+         double cv_hmax_inv;    /* |h| <= 1/hmax_inv                               */\r
+         double cv_etamax;      /* eta <= etamax                                   */\r
+\r
+         /*----------\r
+           Counters \r
+           ----------*/\r
+\r
+         long cv_nst;         /* number of internal steps taken                  */\r
+\r
+         long cv_nfe;         /* number of f calls                               */\r
+         long cv_nfQe;        /* number of fQ calls                              */\r
+         long cv_nfSe;        /* number of fS calls                              */\r
+         long cv_nfeS;        /* number of f calls from sensi DQ                 */\r
+         long cv_nfQSe;       /* number of fQS calls                             */\r
+         long cv_nfQeS;       /* number of fQ calls from sensi DQ                */\r
+\r
+\r
+         long cv_ncfn;        /* number of corrector convergence failures        */\r
+         long cv_ncfnS;       /* number of total sensi. corr. conv. failures     */\r
+         long cv_ncfnS1[];     /* number of sensi. corrector conv. failures       */\r
+\r
+         long cv_nni;         /* number of nonlinear iterations performed        */\r
+         long cv_nniS;        /* number of total sensi. nonlinear iterations     */\r
+         long cv_nniS1[];      /* number of sensi. nonlinear iterations           */\r
+\r
+         long cv_netf;        /* number of error test failures                   */\r
+         long cv_netfQ;       /* number of quadr. error test failures            */\r
+         long cv_netfS;       /* number of sensi. error test failures            */\r
+         long cv_netfQS;      /* number of quadr. sensi. error test failures     */\r
+\r
+         long cv_nsetups;     /* number of setup calls                           */\r
+         long cv_nsetupsS;    /* number of setup calls due to sensitivities      */\r
+\r
+         int cv_nhnil;            /* number of messages issued to the user that\r
+                                     t + h == t for the next iternal step            */\r
+\r
+         /*-----------------------------\r
+           Space requirements for CVODES \r
+           -----------------------------*/\r
+\r
+         long cv_lrw1;        /* no. of realtype 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_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_liw;             /* no. of integer words in CVODES work vectors     */\r
+\r
+         /*----------------\r
+           Step size ratios\r
+           ----------------*/\r
+\r
+         double cv_etaqm1;      /* ratio of new to old h for order q-1             */\r
+         double cv_etaq;        /* ratio of new to old h for order q               */\r
+         double cv_etaqp1;      /* ratio of new to old h for order q+1             */\r
+\r
+         /*------------------\r
+           Linear Solver Data \r
+           ------------------*/\r
+\r
+         /* Linear Solver functions to be called */\r
+\r
+         //int (*cv_linit)(struct CVodeMemRec *cv_mem);\r
+\r
+         //int (*cv_lsetup)(struct CVodeMemRec *cv_mem, int convfail, \r
+                          //N_Vector ypred, N_Vector fpred, booleantype *jcurPtr, \r
+                          //N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3); \r
+\r
+         //int (*cv_lsolve)(struct CVodeMemRec *cv_mem, N_Vector b, N_Vector weight,\r
+                          //N_Vector ycur, N_Vector fcur);\r
+\r
+         //int (*cv_lfree)(struct CVodeMemRec *cv_mem);\r
+\r
+         /* Linear Solver specific memory */\r
+\r
+         //void *cv_lmem;           \r
+\r
+         /* Flag to request a call to the setup routine */\r
+\r
+         boolean cv_forceSetup;\r
+\r
+         /*------------\r
+           Saved Values\r
+           ------------*/\r
+\r
+         int cv_qu;                   /* last successful q value used                */\r
+         long cv_nstlp;           /* step number of last setup call              */\r
+         double cv_h0u;             /* actual initial stepsize                     */\r
+         double cv_hu;              /* last successful h value used                */\r
+         double cv_saved_tq5;       /* saved value of tq[5]                        */\r
+         boolean cv_jcur;         /* is Jacobian info for linear solver current? */\r
+         int cv_convfail;             /* flag storing previous solver failure mode   */\r
+         double cv_tolsf;           /* tolerance scale factor                      */\r
+         int cv_qmax_alloc;           /* qmax used when allocating mem               */\r
+         int cv_qmax_allocQ;          /* qmax used when allocating quad. mem         */\r
+         int cv_qmax_allocS;          /* qmax used when allocating sensi. mem        */\r
+         int cv_qmax_allocQS;         /* qmax used when allocating quad. sensi. mem  */\r
+         int cv_indx_acor;            /* index of zn vector in which acor is saved   */\r
+\r
+         /*--------------------------------------------------------------------\r
+           Flags turned ON by CVodeInit, CVodeSensMalloc, and CVodeQuadMalloc \r
+           and read by CVodeReInit, CVodeSensReInit, and CVodeQuadReInit\r
+           --------------------------------------------------------------------*/\r
+\r
+         boolean cv_VabstolMallocDone;\r
+         boolean cv_MallocDone;\r
+\r
+         boolean cv_VabstolQMallocDone;\r
+         boolean cv_QuadMallocDone;\r
+\r
+         boolean cv_VabstolSMallocDone;\r
+         boolean cv_SabstolSMallocDone;\r
+         boolean cv_SensMallocDone;\r
+\r
+         boolean cv_VabstolQSMallocDone;\r
+         boolean cv_SabstolQSMallocDone;\r
+         boolean cv_QuadSensMallocDone;\r
+\r
+         /*-------------------------------------------\r
+           Error handler function and error ouput file \r
+           -------------------------------------------*/\r
+\r
+         //CVErrHandlerFn cv_ehfun;    /* Error messages are handled by ehfun          */\r
+         CVodeMemRec cv_eh_data;           /* dats pointer passed to ehfun                 */\r
+         RandomAccessFile cv_errfp;             /* CVODES error messages are sent to errfp      */    \r
+\r
+         /*-------------------------\r
+           Stability Limit Detection\r
+           -------------------------*/\r
+\r
+         boolean cv_sldeton;     /* Is Stability Limit Detection on?             */\r
+         double cv_ssdat[][] = new double[6][4];    /* scaled data array for STALD                  */\r
+         int cv_nscon;               /* counter for STALD method                     */\r
+         long cv_nor;            /* counter for number of order reductions       */\r
+\r
+         /*----------------\r
+           Rootfinding Data\r
+           ----------------*/\r
+\r
+         //CVRootFn cv_gfun;        /* Function g for roots sought                     */\r
+         int cv_nrtfn;            /* number of components of g                       */\r
+         int cv_iroots[];          /* array for root information                      */\r
+         int cv_rootdir[];         /* array specifying direction of zero-crossing     */\r
+         double cv_tlo;         /* nearest endpoint of interval in root search     */\r
+         double cv_thi;         /* farthest endpoint of interval in root search    */\r
+         double cv_trout;       /* t value returned by rootfinding routine         */\r
+         double cv_glo[];        /* saved array of g values at t = tlo              */\r
+         double cv_ghi[];        /* saved array of g values at t = thi              */\r
+         double cv_grout[];      /* array of g values at t = trout                  */\r
+         double cv_toutc;       /* copy of tout (if NORMAL mode)                   */\r
+         double cv_ttol;        /* tolerance on root location trout                */\r
+         int cv_taskc;            /* copy of parameter itask                         */\r
+         int cv_irfnd;            /* flag showing whether last step had a root       */\r
+         long cv_nge;         /* counter for g evaluations                       */\r
+         boolean cv_gactive[]; /* array with active/inactive event functions      */\r
+         int cv_mxgnull;          /* number of warning messages about possible g==0  */\r
+\r
+         /*------------------------\r
+           Adjoint sensitivity data\r
+           ------------------------*/\r
+\r
+         boolean cv_adj;             /* SUNTRUE if performing ASA                */\r
+\r
+         CVodeMemRec cv_adj_mem; /* Pointer to adjoint memory structure      */\r
+\r
+         boolean cv_adjMallocDone;\r
+\r
+       } ;\r
+\r
+\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
+       \r
+        // CVodeCreate creates an internal memory block for a problem to\r
+        // be solved by CVODES.\r
+        //\r
+        // lmm  is the type of linear multistep method to be used.\r
+        //      The legal values are CV_ADAMS and CV_BDF (see previous\r
+        //      description).\r
+        //\r
+        // iter  is the type of iteration used to solve the nonlinear\r
+        //       system that arises during each internal time step.\r
+        //       The legal values are CV_FUNCTIONAL and CV_NEWTON.\r
+        //\r
+        //If successful, CVodeCreate returns a pointer to initialized\r
+        // problem memory. This pointer should be passed to CVodeInit.\r
+        // If an initialization error occurs, CVodeCreate prints an error\r
+        // message to standard err and returns NULL.\r
+\r
+     private CVodeMemRec CVodeCreate(int lmm, int iter) {\r
+        int maxord;\r
+         CVodeMemRec cv_mem;\r
+\r
+         /* Test inputs */\r
+\r
+         if ((lmm != CV_ADAMS) && (lmm != CV_BDF)) {\r
+           cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_LMM);\r
+           return null;\r
+         }\r
+         \r
+         if ((iter != CV_FUNCTIONAL) && (iter != CV_NEWTON)) {\r
+           cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_ITER);\r
+           return null;\r
+         }\r
+\r
+         cv_mem = null;\r
+         cv_mem = new CVodeMemRec();\r
+         if (cv_mem == null) {\r
+           cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_CVMEM_FAIL);\r
+           return null;\r
+         }\r
+\r
+         /* Zero out cv_mem */\r
+         //memset(cv_mem, 0, sizeof(struct CVodeMemRec));\r
+\r
+         maxord = (lmm == CV_ADAMS) ? ADAMS_Q_MAX : BDF_Q_MAX;\r
+\r
+         /* copy input parameters into cv_mem */\r
+\r
+         cv_mem.cv_lmm  = lmm;\r
+         cv_mem.cv_iter = iter;\r
+\r
+         /* Set uround */\r
+\r
+         cv_mem.cv_uround = UNIT_ROUNDOFF;\r
+\r
+         /* Set default values for integrator optional inputs */\r
+\r
+         //cv_mem.cv_f          = null;\r
+         //cv_mem.cv_user_data  = null;\r
+         cv_mem.cv_itol       = CV_NN;\r
+         cv_mem.cv_user_efun  = false;\r
+         //cv_mem.cv_efun       = null;\r
+         //cv_mem.cv_e_data     = null;\r
+         //cv_mem.cv_ehfun      = cvErrHandler;\r
+         cv_mem.cv_eh_data    = cv_mem;\r
+         //cv_mem.cv_errfp      = stderr;\r
+         cv_mem.cv_qmax       = maxord;\r
+         cv_mem.cv_mxstep     = MXSTEP_DEFAULT;\r
+         cv_mem.cv_mxhnil     = MXHNIL_DEFAULT;\r
+         cv_mem.cv_sldeton    = false;\r
+         cv_mem.cv_hin        = ZERO;\r
+         cv_mem.cv_hmin       = HMIN_DEFAULT;\r
+         cv_mem.cv_hmax_inv   = HMAX_INV_DEFAULT;\r
+         cv_mem.cv_tstopset   = false;\r
+         cv_mem.cv_maxcor     = NLS_MAXCOR;\r
+         cv_mem.cv_maxnef     = MXNEF;\r
+         cv_mem.cv_maxncf     = MXNCF;\r
+         cv_mem.cv_nlscoef    = CORTES;\r
+\r
+         /* Initialize root finding variables */\r
+\r
+         cv_mem.cv_glo        = null;\r
+         cv_mem.cv_ghi        = null;\r
+         cv_mem.cv_grout      = null;\r
+         cv_mem.cv_iroots     = null;\r
+         cv_mem.cv_rootdir    = null;\r
+         //cv_mem.cv_gfun       = null;\r
+         cv_mem.cv_nrtfn      = 0;  \r
+         cv_mem.cv_gactive    = null;\r
+         cv_mem.cv_mxgnull    = 1;\r
+\r
+         /* Set default values for quad. optional inputs */\r
+\r
+         cv_mem.cv_quadr      = false;\r
+         //cv_mem.cv_fQ         = null;\r
+         cv_mem.cv_errconQ    = false;\r
+         cv_mem.cv_itolQ      = CV_NN;\r
+\r
+         /* 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         = cvSensRhsInternalDQ;\r
+         //cv_mem.cv_fS1        = cvSensRhs1InternalDQ;\r
+         cv_mem.cv_fSDQ       = true;\r
+         cv_mem.cv_ifS        = CV_ONESENS;\r
+         cv_mem.cv_DQtype     = CV_CENTERED;\r
+         cv_mem.cv_DQrhomax   = ZERO;\r
+         cv_mem.cv_p          = null;\r
+         cv_mem.cv_pbar       = null;\r
+         cv_mem.cv_plist      = null;\r
+         cv_mem.cv_errconS    = false;\r
+         cv_mem.cv_maxcorS    = NLS_MAXCOR;\r
+         cv_mem.cv_ncfS1      = null;\r
+         cv_mem.cv_ncfnS1     = null;\r
+         cv_mem.cv_nniS1      = null;\r
+         cv_mem.cv_itolS      = CV_NN;\r
+\r
+         /* Set default values for quad. sensi. optional inputs */\r
+\r
+         cv_mem.cv_quadr_sensi = false;\r
+         //cv_mem.cv_fQS         = null;\r
+         //cv_mem.cv_fQS_data    = null;\r
+         cv_mem.cv_fQSDQ       = true;\r
+         cv_mem.cv_errconQS    = false;\r
+         cv_mem.cv_itolQS      = CV_NN;\r
+\r
+         /* Set default for ASA */\r
+\r
+         cv_mem.cv_adj         = false;\r
+         cv_mem.cv_adj_mem     = null;\r
+\r
+         /* Set the saved values for qmax_alloc */\r
+\r
+         cv_mem.cv_qmax_alloc  = maxord;\r
+         cv_mem.cv_qmax_allocQ = maxord;\r
+         cv_mem.cv_qmax_allocS = maxord;\r
+\r
+         /* Initialize lrw and liw */\r
+\r
+         cv_mem.cv_lrw = 65 + 2*L_MAX + NUM_TESTS;\r
+         cv_mem.cv_liw = 52;\r
+\r
+         /* No mallocs have been done yet */\r
+\r
+         cv_mem.cv_VabstolMallocDone   = false;\r
+         cv_mem.cv_MallocDone          = false;\r
+\r
+         cv_mem.cv_VabstolQMallocDone  = false;\r
+         cv_mem.cv_QuadMallocDone      = false;\r
+\r
+         cv_mem.cv_VabstolSMallocDone  = false;\r
+         cv_mem.cv_SabstolSMallocDone  = false;\r
+         cv_mem.cv_SensMallocDone      = false;\r
+\r
+         cv_mem.cv_VabstolQSMallocDone = false;\r
+         cv_mem.cv_SabstolQSMallocDone = false;\r
+         cv_mem.cv_QuadSensMallocDone  = false;\r
+\r
+         cv_mem.cv_adjMallocDone       = false;\r
+\r
+         /* Return pointer to CVODES memory block */\r
+\r
+         return cv_mem;\r
+\r
+     }\r
+     \r
+     /*\r
+      * cvProcessError is a high level error handling function.\r
+      * - If cv_mem==NULL it prints the error message to stderr.\r
+      * - Otherwise, it sets up and calls the error handling function \r
+      *   pointed to by cv_ehfun.\r
+      */\r
+\r
+     void cvProcessError(CVodeMemRec cv_mem, \r
+                         int error_code,  String module, String fname, \r
+                         String msgfmt)\r
+     {\r
+        MipavUtil.displayError("In module " + module + " in function " + fname + " msgfmt");\r
+       //va_list ap;\r
+       //char msg[256];\r
+\r
+       /* Initialize the argument pointer variable \r
+          (msgfmt is the last required argument to cvProcessError) */\r
+\r
+       //va_start(ap, msgfmt);\r
+\r
+       /* Compose the message */\r
+\r
+       //vsprintf(msg, msgfmt, ap);\r
+\r
+       //if (cv_mem == NULL) {    /* We write to stderr */\r
+     //#ifndef NO_FPRINTF_OUTPUT\r
+         //fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);\r
+         //fprintf(stderr, "%s\n\n", msg);\r
+     //#endif\r
+\r
+       //} else {                 /* We can call ehfun */\r
+         //cv_mem->cv_ehfun(error_code, module, fname, msg, cv_mem->cv_eh_data);\r
+       //}\r
+\r
+       /* Finalize argument processing */\r
+       //va_end(ap);\r
+\r
+       return;\r
+     }\r
+     \r
+     /*-----------------------------------------------------------------*/\r
+\r
+     /*\r
+      * CVodeInit\r
+      * \r
+      * CVodeInit allocates and initializes memory for a problem. All \r
+      * problem inputs are checked for errors. If any error occurs during \r
+      * initialization, it is reported to the file whose file pointer is \r
+      * errfp and an error flag is returned. Otherwise, it returns CV_SUCCESS\r
+      */\r
+\r
+     int CVodeInit(CVodeMemRec cv_mem, double t0, double y0[])\r
+     {\r
+       boolean nvectorOK, allocOK;\r
+       long lrw1, liw1;\r
+       int i,k;\r
+\r
+       /* Check cvode_mem */\r
+\r
+       if (cv_mem==null) {\r
+         cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeInit",\r
+                        MSGCV_NO_MEM);\r
+         return(CV_MEM_NULL);\r
+       }\r
+\r
+       /* Check for legal input parameters */\r
+\r
+       if (y0==null) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
+                        MSGCV_NULL_Y0);\r
+         return(CV_ILL_INPUT);\r
+       }\r
+\r
+       //if (f == NULL) {\r
+         //cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
+                        //MSGCV_NULL_F);\r
+         //return(CV_ILL_INPUT);\r
+       //}\r
+\r
+       /* Test if all required vector operations are implemented */\r
+\r
+       /*nvectorOK = cvCheckNvector(y0);\r
+       if(!nvectorOK) {\r
+         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
+                        MSGCV_BAD_NVECTOR);\r
+         return(CV_ILL_INPUT);\r
+       }*/\r
+\r
+       /* Set space requirements for one N_Vector */\r
+\r
+       /*if (y0->ops->nvspace != NULL) {\r
+         N_VSpace(y0, &lrw1, &liw1);\r
+       } else {\r
+         lrw1 = 0;\r
+         liw1 = 0;\r
+       }\r
+       cv_mem->cv_lrw1 = lrw1;\r
+       cv_mem->cv_liw1 = liw1;*/\r
+\r
+       /* Allocate the vectors (using y0 as a template) */\r
+\r
+      /* allocOK = cvAllocVectors(cv_mem, y0);\r
+       if (!allocOK) {\r
+         cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeInit",\r
+                        MSGCV_MEM_FAIL);\r
+         return(CV_MEM_FAIL);\r
+       }*/\r
+\r
+       /* All error checking is complete at this point */\r
+\r
+       /* Copy the input parameters into CVODES state */\r
+\r
+       /*cv_mem->cv_f  = f;\r
+       cv_mem->cv_tn = t0;*/\r
+\r
+       /* Set step parameters */\r
+\r
+       /*cv_mem->cv_q      = 1;\r
+       cv_mem->cv_L      = 2;\r
+       cv_mem->cv_qwait  = cv_mem->cv_L;\r
+       cv_mem->cv_etamax = ETAMX1;\r
+\r
+       cv_mem->cv_qu     = 0;\r
+       cv_mem->cv_hu     = ZERO;\r
+       cv_mem->cv_tolsf  = ONE;*/\r
+\r
+       /* Set the linear solver addresses to NULL.\r
+          (We check != NULL later, in CVode, if using CV_NEWTON.) */\r
+\r
+       /*cv_mem->cv_linit  = NULL;\r
+       cv_mem->cv_lsetup = NULL;\r
+       cv_mem->cv_lsolve = NULL;\r
+       cv_mem->cv_lfree  = NULL;\r
+       cv_mem->cv_lmem   = NULL;*/\r
+\r
+       /* Set forceSetup to SUNFALSE */\r
+\r
+       //cv_mem->cv_forceSetup = SUNFALSE;\r
+\r
+       /* Initialize zn[0] in the history array */\r
+\r
+       //N_VScale(ONE, y0, cv_mem->cv_zn[0]);\r
+\r
+       /* Initialize all the counters */\r
+\r
+       /*cv_mem->cv_nst     = 0;\r
+       cv_mem->cv_nfe     = 0;\r
+       cv_mem->cv_ncfn    = 0;\r
+       cv_mem->cv_netf    = 0;\r
+       cv_mem->cv_nni     = 0;\r
+       cv_mem->cv_nsetups = 0;\r
+       cv_mem->cv_nhnil   = 0;\r
+       cv_mem->cv_nstlp   = 0;\r
+       cv_mem->cv_nscon   = 0;\r
+       cv_mem->cv_nge     = 0;\r
+\r
+       cv_mem->cv_irfnd   = 0;*/\r
+\r
+       /* Initialize other integrator optional outputs */\r
+\r
+       /*cv_mem->cv_h0u      = ZERO;\r
+       cv_mem->cv_next_h   = ZERO;\r
+       cv_mem->cv_next_q   = 0;*/\r
+\r
+       /* Initialize Stablilty Limit Detection data */\r
+       /* NOTE: We do this even if stab lim det was not\r
+          turned on yet. This way, the user can turn it\r
+          on at any time */\r
+\r
+       /*cv_mem->cv_nor = 0;\r
+       for (i = 1; i <= 5; i++)\r
+         for (k = 1; k <= 3; k++) \r
+           cv_mem->cv_ssdat[i-1][k-1] = ZERO;*/\r
+\r
+       /* Problem has been successfully initialized */\r
+\r
+       /*cv_mem->cv_MallocDone = SUNTRUE;*/\r
+\r
+       return(CV_SUCCESS);\r
+     }\r
+     \r
+     /*\r
+      * cvCheckNvector\r
+      * This routine checks if all required vector operations are present.\r
+      * If any of them is missing it returns SUNFALSE.\r
+      */\r
+\r
+     /*static booleantype cvCheckNvector(N_Vector tmpl)\r
+     {\r
+       if((tmpl->ops->nvclone     == NULL) || // clone\r
+          (tmpl->ops->nvdestroy   == NULL) || // set null\r
+          (tmpl->ops->nvlinearsum == NULL) || // (double a, double x[], double b, double y[], double z[])\r
+                                              // z[i] = a*x[i] + b*y[i];\r
+          (tmpl->ops->nvconst     == NULL) ||\r
+          (tmpl->ops->nvprod      == NULL) ||\r
+          (tmpl->ops->nvdiv       == NULL) ||\r
+          (tmpl->ops->nvscale     == NULL) ||\r
+          (tmpl->ops->nvabs       == NULL) ||\r
+          (tmpl->ops->nvinv       == NULL) ||\r
+          (tmpl->ops->nvaddconst  == NULL) ||\r
+          (tmpl->ops->nvmaxnorm   == NULL) ||\r
+          (tmpl->ops->nvwrmsnorm  == NULL) ||\r
+          (tmpl->ops->nvmin       == NULL))\r
+         return(SUNFALSE);\r
+       else\r
+         return(SUNTRUE);\r
+     }*/\r
+\r
+\r
+\r
+}
\ No newline at end of file