8b08cad87ac024a02f197baf22ea9bd70e66a637
[mipav.git] / mipav / src / gov / nih / mipav / model / algorithms / CVODES.java
1 package gov.nih.mipav.model.algorithms;\r
2 \r
3 import java.io.RandomAccessFile;\r
4 \r
5 import gov.nih.mipav.view.MipavUtil;\r
6 import gov.nih.mipav.view.Preferences;\r
7 \r
8 public abstract class CVODES {\r
9         // This is a port of code from CVODES, a stiff and non-stiff ODE solver, from C to Java\r
10         /*Copyright (c) 2002-2016, Lawrence Livermore National Security. \r
11         Produced at the Lawrence Livermore National Laboratory.\r
12         Written by A.C. Hindmarsh, D.R. Reynolds, R. Serban, C.S. Woodward,\r
13         S.D. Cohen, A.G. Taylor, S. Peles, L.E. Banks, and D. Shumaker.\r
14         LLNL-CODE-667205    (ARKODE)\r
15         UCRL-CODE-155951    (CVODE)\r
16         UCRL-CODE-155950    (CVODES)\r
17         UCRL-CODE-155952    (IDA)\r
18         UCRL-CODE-237203    (IDAS)\r
19         LLNL-CODE-665877    (KINSOL)\r
20         All rights reserved. \r
21 \r
22         This file is part of SUNDIALS.  For details, \r
23         see http://computation.llnl.gov/projects/sundials\r
24 \r
25         Redistribution and use in source and binary forms, with or without\r
26         modification, are permitted provided that the following conditions\r
27         are met:\r
28 \r
29         1. Redistributions of source code must retain the above copyright\r
30         notice, this list of conditions and the disclaimer below.\r
31 \r
32         2. Redistributions in binary form must reproduce the above copyright\r
33         notice, this list of conditions and the disclaimer (as noted below)\r
34         in the documentation and/or other materials provided with the\r
35         distribution.\r
36 \r
37         3. Neither the name of the LLNS/LLNL nor the names of its contributors\r
38         may be used to endorse or promote products derived from this software\r
39         without specific prior written permission.\r
40 \r
41         THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \r
42         "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT \r
43         LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS \r
44         FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL \r
45         LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF \r
46         ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, \r
47         SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED \r
48         TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, \r
49         DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY \r
50         THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
51         (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
52         OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
53 \r
54         Additional BSD Notice\r
55         ---------------------\r
56         1. This notice is required to be provided under our contract with\r
57         the U.S. Department of Energy (DOE). This work was produced at\r
58         Lawrence Livermore National Laboratory under Contract \r
59         No. DE-AC52-07NA27344 with the DOE.\r
60 \r
61         2. Neither the United States Government nor Lawrence Livermore \r
62         National Security, LLC nor any of their employees, makes any warranty, \r
63         express or implied, or assumes any liability or responsibility for the\r
64         accuracy, completeness, or usefulness of any */\r
65         \r
66         /* Basic CVODES constants */\r
67 \r
68         final int ADAMS_Q_MAX = 12;      /* max value of q for lmm == ADAMS    */\r
69         final int BDF_Q_MAX = 5;      /* max value of q for lmm == BDF      */\r
70         final int  Q_MAX = ADAMS_Q_MAX;  /* max value of q for either lmm      */\r
71         final int L_MAX  = (Q_MAX+1);    /* max value of L for either lmm      */\r
72         final int NUM_TESTS = 5;      /* number of error test quantities    */\r
73         \r
74     double DBL_EPSILON;\r
75     double UNIT_ROUNDOFF;\r
76 \r
77 \r
78 \r
79          // lmm:   The user of the CVODES package specifies whether to use\r
80          // the CV_ADAMS or CV_BDF (backward differentiation formula)\r
81          // linear multistep method. The BDF method is recommended\r
82          // for stiff problems, and the CV_ADAMS method is recommended\r
83          // for nonstiff problems.\r
84         final int CV_ADAMS = 1;\r
85         final int CV_BDF = 2;\r
86         \r
87         //iter:  At each internal time step, a nonlinear equation must\r
88         //        be solved. The user can specify either CV_FUNCTIONAL\r
89         //        iteration, which does not require linear algebra, or a\r
90         //        CV_NEWTON iteration, which requires the solution of linear\r
91         //        systems. In the CV_NEWTON case, the user also specifies a\r
92         //        CVODE linear solver. CV_NEWTON is recommended in case of\r
93         //        stiff problems.\r
94         final int CV_FUNCTIONAL = 1;\r
95         final int CV_NEWTON = 2;\r
96         \r
97         /*\r
98          * Control constants for type of sensitivity RHS\r
99          * ---------------------------------------------\r
100          */\r
101 \r
102         final int CV_ONESENS = 1;\r
103         final int CV_ALLSENS = 2;\r
104 \r
105         \r
106         /*\r
107          * Control constants for tolerances\r
108          * --------------------------------\r
109          */\r
110 \r
111         final int CV_NN = 0;\r
112         final int CV_SS = 1;\r
113         final int CV_SV = 2;\r
114         final int CV_WF = 3;\r
115         final int CV_EE = 4;\r
116         \r
117         /* DQtype */\r
118         final int CV_CENTERED = 1;\r
119         final int CV_FORWARD = 2;\r
120         \r
121         /* \r
122          * ----------------------------------------\r
123          * CVODES return flags\r
124          * ----------------------------------------\r
125          */\r
126 \r
127         final int  CV_SUCCESS = 0;\r
128         final int CV_TSTOP_RETURN = 1;\r
129         final int CV_ROOT_RETURN = 2;\r
130 \r
131         final int CV_WARNING = 99;\r
132 \r
133         final int CV_TOO_MUCH_WORK = -1;\r
134         final int CV_TOO_MUCH_ACC = -2;\r
135         final int CV_ERR_FAILURE = -3;\r
136         final int CV_CONV_FAILURE = -4;\r
137 \r
138         final int CV_LINIT_FAIL = -5;\r
139         final int CV_LSETUP_FAIL = -6;\r
140         final int CV_LSOLVE_FAIL = -7;\r
141         final int CV_RHSFUNC_FAIL = -8;\r
142         final int CV_FIRST_RHSFUNC_ERR = -9;\r
143         final int CV_REPTD_RHSFUNC_ERR = -10;\r
144         final int CV_UNREC_RHSFUNC_ERR = -11;\r
145         final int  CV_RTFUNC_FAIL = -12;\r
146 \r
147         final int CV_MEM_FAIL = -20;\r
148         final int CV_MEM_NULL = -21;\r
149         final int CV_ILL_INPUT = -22;\r
150         final int CV_NO_MALLOC = -23;\r
151         final int CV_BAD_K = -24;\r
152         final int CV_BAD_T = -25;\r
153         final int CV_BAD_DKY = -26;\r
154         final int CV_TOO_CLOSE = -27;\r
155 \r
156         final int CV_NO_QUAD = -30;\r
157         final int CV_QRHSFUNC_FAIL = -31;\r
158         final int CV_FIRST_QRHSFUNC_ERR = -32;\r
159         final int CV_REPTD_QRHSFUNC_ERR = -33;\r
160         final int CV_UNREC_QRHSFUNC_ERR = -34;\r
161 \r
162         final int CV_NO_SENS = -40;\r
163         final int CV_SRHSFUNC_FAIL = -41;\r
164         final int CV_FIRST_SRHSFUNC_ERR = -42;\r
165         final int CV_REPTD_SRHSFUNC_ERR = -43;\r
166         final int CV_UNREC_SRHSFUNC_ERR = -44;\r
167 \r
168         final int CV_BAD_IS = -45;\r
169 \r
170         final int CV_NO_QUADSENS = -50;\r
171         final int CV_QSRHSFUNC_FAIL = -51;\r
172         final int CV_FIRST_QSRHSFUNC_ERR = -52;\r
173         final int CV_REPTD_QSRHSFUNC_ERR = -53;\r
174         final int CV_UNREC_QSRHSFUNC_ERR = -54;\r
175 \r
176         final double HMIN_DEFAULT = 0.0;    /* hmin default value     */\r
177         final double HMAX_INV_DEFAULT = 0.0;    /* hmax_inv default value */\r
178         final int MXHNIL_DEFAULT =  10;             /* mxhnil default value   */\r
179         final long MXSTEP_DEFAULT = 500L;            /* mxstep default value   */\r
180         final int NLS_MAXCOR = 3;\r
181         final int MXNCF = 10;\r
182         final int MXNEF = 7;\r
183         final double CORTES = 0.1;\r
184         final double ZERO = 0.0;\r
185         final double ONE = 1.0;\r
186         final double ETAMX1 = 10000.0; \r
187 \r
188         \r
189         final String MSG_TIME = "t = %lg";\r
190         final String MSG_TIME_H = "t = %lg and h = %lg";\r
191         final String MSG_TIME_INT = "t = %lg is not between tcur - hu = %lg and tcur = %lg.";\r
192         final String MSG_TIME_TOUT = "tout = %lg";\r
193         final String MSG_TIME_TSTOP = "tstop = %lg";\r
194 \r
195         \r
196         /* Initialization and I/O error messages */\r
197 \r
198         final String MSGCV_NO_MEM = "cvode_mem = NULL illegal.";\r
199         final String MSGCV_CVMEM_FAIL = "Allocation of cvode_mem failed.";\r
200         final String MSGCV_MEM_FAIL = "A memory request failed.";\r
201         final String MSGCV_BAD_LMM = "Illegal value for lmm. The legal values are CV_ADAMS and CV_BDF.";\r
202         final String MSGCV_BAD_ITER = "Illegal value for iter. The legal values are CV_FUNCTIONAL and CV_NEWTON.";\r
203         final String MSGCV_NO_MALLOC = "Attempt to call before CVodeInit.";\r
204         final String MSGCV_NEG_MAXORD = "maxord <= 0 illegal.";\r
205         final String MSGCV_BAD_MAXORD = "Illegal attempt to increase maximum method order.";\r
206         final String MSGCV_SET_SLDET = "Attempt to use stability limit detection with the CV_ADAMS method illegal.";\r
207         final String MSGCV_NEG_HMIN = "hmin < 0 illegal.";\r
208         final String MSGCV_NEG_HMAX = "hmax < 0 illegal.";\r
209         final String MSGCV_BAD_HMIN_HMAX = "Inconsistent step size limits: hmin > hmax.";\r
210         final String MSGCV_BAD_RELTOL = "reltol < 0 illegal.";\r
211         final String MSGCV_BAD_ABSTOL = "abstol has negative component(s) (illegal).";\r
212         final String MSGCV_NULL_ABSTOL = "abstol = NULL illegal.";\r
213         final String MSGCV_NULL_Y0 = "y0 = NULL illegal.";\r
214         final String MSGCV_NULL_F = "f = NULL illegal.";\r
215         final String MSGCV_NULL_G = "g = NULL illegal.";\r
216         final String MSGCV_BAD_NVECTOR = "A required vector operation is not implemented.";\r
217         final String MSGCV_BAD_K = "Illegal value for k.";\r
218         final String MSGCV_NULL_DKY = "dky = NULL illegal.";\r
219         final String MSGCV_BAD_T = "Illegal value for t.";\r
220         final String MSGCV_NO_ROOT = "Rootfinding was not initialized.";\r
221 \r
222         final String MSGCV_NO_QUAD  = "Quadrature integration not activated.";\r
223         final String MSGCV_BAD_ITOLQ = "Illegal value for itolQ. The legal values are CV_SS and CV_SV.";\r
224         final String MSGCV_NULL_ABSTOLQ = "abstolQ = NULL illegal.";\r
225         final String MSGCV_BAD_RELTOLQ = "reltolQ < 0 illegal.";\r
226         final String MSGCV_BAD_ABSTOLQ = "abstolQ has negative component(s) (illegal).";  \r
227 \r
228         final String MSGCV_SENSINIT_2 = "Sensitivity analysis already initialized.";\r
229         final String MSGCV_NO_SENSI  = "Forward sensitivity analysis not activated.";\r
230         final String MSGCV_BAD_ITOLS = "Illegal value for itolS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
231         final String MSGCV_NULL_ABSTOLS = "abstolS = NULL illegal.";\r
232         final String MSGCV_BAD_RELTOLS = "reltolS < 0 illegal.";\r
233         final String MSGCV_BAD_ABSTOLS = "abstolS has negative component(s) (illegal).";  \r
234         final String MSGCV_BAD_PBAR = "pbar has zero component(s) (illegal).";\r
235         final String MSGCV_BAD_PLIST = "plist has negative component(s) (illegal).";\r
236         final String MSGCV_BAD_NS = "NS <= 0 illegal.";\r
237         final String MSGCV_NULL_YS0 = "yS0 = NULL illegal.";\r
238         final String MSGCV_BAD_ISM = "Illegal value for ism. Legal values are: CV_SIMULTANEOUS, CV_STAGGERED and CV_STAGGERED1.";\r
239         final String MSGCV_BAD_IFS = "Illegal value for ifS. Legal values are: CV_ALLSENS and CV_ONESENS.";\r
240         final String MSGCV_BAD_ISM_IFS = "Illegal ism = CV_STAGGERED1 for CVodeSensInit.";\r
241         final String MSGCV_BAD_IS = "Illegal value for is.";\r
242         final String MSGCV_NULL_DKYA = "dkyA = NULL illegal.";\r
243         final String MSGCV_BAD_DQTYPE = "Illegal value for DQtype. Legal values are: CV_CENTERED and CV_FORWARD.";\r
244         final String MSGCV_BAD_DQRHO = "DQrhomax < 0 illegal.";\r
245 \r
246         final String MSGCV_BAD_ITOLQS = "Illegal value for itolQS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
247         final String MSGCV_NULL_ABSTOLQS = "abstolQS = NULL illegal.";\r
248         final String MSGCV_BAD_RELTOLQS = "reltolQS < 0 illegal.";\r
249         final String MSGCV_BAD_ABSTOLQS = "abstolQS has negative component(s) (illegal).";  \r
250         final String MSGCV_NO_QUADSENSI = "Forward sensitivity analysis for quadrature variables not activated.";\r
251         final String MSGCV_NULL_YQS0 = "yQS0 = NULL illegal.";\r
252 \r
253         /* CVode Error Messages */\r
254 \r
255         final String MSGCV_NO_TOL = "No integration tolerances have been specified.";\r
256         final String MSGCV_LSOLVE_NULL = "The linear solver's solve routine is NULL.";\r
257         final String MSGCV_YOUT_NULL = "yout = NULL illegal.";\r
258         final String MSGCV_TRET_NULL = "tret = NULL illegal.";\r
259         final String MSGCV_BAD_EWT = "Initial ewt has component(s) equal to zero (illegal).";\r
260         final String MSGCV_EWT_NOW_BAD = "At " + MSG_TIME + ", a component of ewt has become <= 0.";\r
261         final String MSGCV_BAD_ITASK = "Illegal value for itask.";\r
262         final String MSGCV_BAD_H0 = "h0 and tout - t0 inconsistent.";\r
263         final String MSGCV_BAD_TOUT = "Trouble interpolating at " + MSG_TIME_TOUT + ". tout too far back in direction of integration";\r
264         final String MSGCV_EWT_FAIL = "The user-provide EwtSet function failed.";\r
265         final String MSGCV_EWT_NOW_FAIL = "At " + MSG_TIME + ", the user-provide EwtSet function failed.";\r
266         final String MSGCV_LINIT_FAIL = "The linear solver's init routine failed.";\r
267         final String MSGCV_HNIL_DONE = "The above warning has been issued mxhnil times and will not be issued again for this problem.";\r
268         final String MSGCV_TOO_CLOSE = "tout too close to t0 to start integration.";\r
269         final String MSGCV_MAX_STEPS = "At " + MSG_TIME + ", mxstep steps taken before reaching tout.";\r
270         final String MSGCV_TOO_MUCH_ACC = "At " + MSG_TIME + ", too much accuracy requested.";\r
271         final String MSGCV_HNIL = "Internal " + MSG_TIME_H + " are such that t + h = t on the next step. The solver will continue anyway.";\r
272         final String MSGCV_ERR_FAILS = "At " + MSG_TIME_H + ", the error test failed repeatedly or with |h| = hmin.";\r
273         final String MSGCV_CONV_FAILS = "At " + MSG_TIME_H + ", the corrector convergence test failed repeatedly or with |h| = hmin.";\r
274         final String MSGCV_SETUP_FAILED = "At " + MSG_TIME + ", the setup routine failed in an unrecoverable manner.";\r
275         final String MSGCV_SOLVE_FAILED = "At " + MSG_TIME + ", the solve routine failed in an unrecoverable manner.";\r
276         final String MSGCV_RHSFUNC_FAILED = "At " + MSG_TIME + ", the right-hand side routine failed in an unrecoverable manner.";\r
277         final String MSGCV_RHSFUNC_UNREC = "At " + MSG_TIME + ", the right-hand side failed in a recoverable manner, but no recovery is possible.";\r
278         final String MSGCV_RHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable right-hand side function errors.";\r
279         final String MSGCV_RHSFUNC_FIRST = "The right-hand side routine failed at the first call.";\r
280         final String MSGCV_RTFUNC_FAILED = "At " + MSG_TIME + ", the rootfinding routine failed in an unrecoverable manner.";\r
281         final String MSGCV_CLOSE_ROOTS = "Root found at and very near " + MSG_TIME + ".";\r
282         final String MSGCV_BAD_TSTOP = "The value " + MSG_TIME_TSTOP + " is behind current " + MSG_TIME + " in the direction of integration.";\r
283         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
284 \r
285         final String MSGCV_NO_TOLQ = "No integration tolerances for quadrature variables have been specified.";\r
286         final String MSGCV_BAD_EWTQ = "Initial ewtQ has component(s) equal to zero (illegal).";\r
287         final String MSGCV_EWTQ_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQ has become <= 0.";\r
288         final String MSGCV_QRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature right-hand side routine failed in an unrecoverable manner.";\r
289         final String MSGCV_QRHSFUNC_UNREC = "At " + MSG_TIME + ", the quadrature right-hand side failed in a recoverable manner, but no recovery is possible.";\r
290         final String MSGCV_QRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature right-hand side function errors.";\r
291         final String MSGCV_QRHSFUNC_FIRST = "The quadrature right-hand side routine failed at the first call.";\r
292 \r
293         final String MSGCV_NO_TOLS = "No integration tolerances for sensitivity variables have been specified.";\r
294         final String MSGCV_NULL_P = "p = NULL when using internal DQ for sensitivity RHS illegal.";\r
295         final String MSGCV_BAD_EWTS = "Initial ewtS has component(s) equal to zero (illegal).";\r
296         final String MSGCV_EWTS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtS has become <= 0.";\r
297         final String MSGCV_SRHSFUNC_FAILED = "At " + MSG_TIME + ", the sensitivity right-hand side routine failed in an unrecoverable manner.";\r
298         final String MSGCV_SRHSFUNC_UNREC = "At " + MSG_TIME + ", the sensitivity right-hand side failed in a recoverable manner, but no recovery is possible.";\r
299         final String MSGCV_SRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable sensitivity right-hand side function errors.";\r
300         final String MSGCV_SRHSFUNC_FIRST = "The sensitivity right-hand side routine failed at the first call.";\r
301 \r
302         final String MSGCV_NULL_FQ = "CVODES is expected to use DQ to evaluate the RHS of quad. sensi., but quadratures were not initialized.";\r
303         final String MSGCV_NO_TOLQS = "No integration tolerances for quadrature sensitivity variables have been specified.";\r
304         final String MSGCV_BAD_EWTQS = "Initial ewtQS has component(s) equal to zero (illegal).";\r
305         final String MSGCV_EWTQS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQS has become <= 0.";\r
306         final String MSGCV_QSRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature sensitivity right-hand side routine failed in an unrecoverable manner.";\r
307         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
308         final String MSGCV_QSRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature sensitivity right-hand side function errors.";\r
309         final String MSGCV_QSRHSFUNC_FIRST = "The quadrature sensitivity right-hand side routine failed at the first call.";\r
310 \r
311         /* \r
312          * =================================================================\r
313          *   C V O D E A    E R R O R    M E S S A G E S\r
314          * =================================================================\r
315          */\r
316 \r
317         final String MSGCV_NO_ADJ = "Illegal attempt to call before calling CVodeAdjMalloc.";\r
318         final String MSGCV_BAD_STEPS = "Steps nonpositive illegal.";\r
319         final String MSGCV_BAD_INTERP = "Illegal value for interp.";\r
320         final String MSGCV_BAD_WHICH  = "Illegal value for which.";\r
321         final String MSGCV_NO_BCK = "No backward problems have been defined yet.";\r
322         final String MSGCV_NO_FWD = "Illegal attempt to call before calling CVodeF.";\r
323         final String MSGCV_BAD_TB0 = "The initial time tB0 for problem %d is outside the interval over which the forward problem was solved.";\r
324         final String MSGCV_BAD_SENSI = "At least one backward problem requires sensitivities, but they were not stored for interpolation.";\r
325         final String MSGCV_BAD_ITASKB = "Illegal value for itaskB. Legal values are CV_NORMAL and CV_ONE_STEP.";\r
326         final String MSGCV_BAD_TBOUT  = "The final time tBout is outside the interval over which the forward problem was solved.";\r
327         final String MSGCV_BACK_ERROR  = "Error occured while integrating backward problem # %d"; \r
328         final String MSGCV_BAD_TINTERP = "Bad t = %g for interpolation.";\r
329         final String MSGCV_WRONG_INTERP = "This function cannot be called for the specified interp type.";\r
330 \r
331 \r
332     final int cvsRoberts_dns = 1;\r
333     int problem = cvsRoberts_dns;\r
334         \r
335         \r
336         public CVODES() {\r
337                 // eps returns the distance from 1.0 to the next larger double-precision\r
338                 // number, that is, eps = 2^-52.\r
339                 // epsilon = D1MACH(4)\r
340                 // Machine epsilon is the smallest positive epsilon such that\r
341                 // (1.0 + epsilon) != 1.0.\r
342                 // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
343                 // epsilon = 2.2204460e-16\r
344                 // epsilon is called the largest relative spacing\r
345                 DBL_EPSILON = 1.0;\r
346                 double neweps = 1.0;\r
347 \r
348                 while (true) {\r
349 \r
350                         if (1.0 == (1.0 + neweps)) {\r
351                                 break;\r
352                         } else {\r
353                                 DBL_EPSILON = neweps;\r
354                                 neweps = neweps / 2.0;\r
355                         }\r
356                 } // while(true)\r
357                 \r
358             UNIT_ROUNDOFF = DBL_EPSILON;\r
359             if (problem == cvsRoberts_dns) {\r
360                 runcvsRoberts_dns();\r
361             }\r
362             return;\r
363         }\r
364         \r
365         private void runcvsRoberts_dns() {\r
366                 /* -----------------------------------------------------------------\r
367                  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and\r
368                  *                Radu Serban @ LLNL\r
369                  * -----------------------------------------------------------------\r
370                  * Example problem:\r
371                  * \r
372                  * The following is a simple example problem, with the coding\r
373                  * needed for its solution by CVODE. The problem is from\r
374                  * chemical kinetics, and consists of the following three rate\r
375                  * equations:         \r
376                  *    dy1/dt = -.04*y1 + 1.e4*y2*y3\r
377                  *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2\r
378                  *    dy3/dt = 3.e7*(y2)^2\r
379                  * on the interval from t = 0.0 to t = 4.e10, with initial\r
380                  * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.\r
381                  * While integrating the system, we also use the rootfinding\r
382                  * feature to find the points at which y1 = 1e-4 or at which\r
383                  * y3 = 0.01. This program solves the problem with the BDF method,\r
384                  * Newton iteration with the SUNDENSE dense linear solver, and a\r
385                  * user-supplied Jacobian routine.\r
386                  * It uses a scalar relative tolerance and a vector absolute\r
387                  * tolerance. Output is printed in decades from t = .4 to t = 4.e10.\r
388                  * Run statistics (optional outputs) are printed at the end.\r
389                  * -----------------------------------------------------------------*/\r
390                 \r
391                 /** Problem Constants */\r
392                 final int NEQ = 3; // Number of equations\r
393                 final double Y1 = 1.0; // Initial y components\r
394                 final double Y2 = 0.0;\r
395                 final double Y3 = 0.0;\r
396                 final double RTOL = 1.0E-4; // scalar relative tolerance\r
397                 final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
398                 final double ATOL2 = 1.0E-14;\r
399                 final double ATOL3 = 1.0E-6;\r
400                 final double T0 = 0.0; // initial time\r
401                 final double T1 = 0.4; // first output time\r
402                 final double TMULT = 10.0; // output time factor\r
403                 final int NOUT = 12; // number of output times\r
404                 int f = cvsRoberts_dns;\r
405                 int g = cvsRoberts_dns;\r
406                 \r
407                 // Set initial conditions\r
408                 NVector y = new NVector();\r
409                 NVector abstol = new NVector();\r
410                 double yr[] = new double[]{Y1,Y2,Y3};\r
411                 N_VNew_Serial(y, NEQ);\r
412                 y.setData(yr);\r
413                 double reltol = RTOL; // Set the scalar relative tolerance\r
414                 // Set the vector absolute tolerance\r
415                 double abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
416                 N_VNew_Serial(abstol,NEQ);\r
417                 abstol.setData(abstolr);\r
418                 CVodeMemRec cvode_mem;\r
419                 int flag;\r
420                 double A[][];\r
421                 SUNLinearSolver LS;\r
422                 \r
423                 // Call CVodeCreate to create the solver memory and specify the\r
424                 // Backward Differentiation Formula and the use of a Newton\r
425                 // iteration\r
426                 cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);\r
427                 if (cvode_mem == null) {\r
428                     return;     \r
429                 }\r
430                 \r
431                 // Call CVodeInit to initialize the integrator memory and specify the\r
432                 // user's right hand side function in y' = f(t,y), the initial time T0, and\r
433             // the initial dependent variable vector y\r
434                 flag = CVodeInit(cvode_mem, f, T0, y);\r
435                 if (flag != CV_SUCCESS) {\r
436                         return;\r
437                 }\r
438                 \r
439                 // Call CVodeSVtolerances to specify the scalar relative tolerance\r
440                 // and vector absolute tolerances\r
441                 flag = CVodeSVtolerances(cvode_mem, reltol, abstol);\r
442                 if (flag != CV_SUCCESS) {\r
443                         return;\r
444                 }\r
445                 \r
446                 // Call CVRootInit to specify the root function g with 2 components\r
447                 flag = CVodeRootInit(cvode_mem, 2, g);\r
448                 if (flag != CV_SUCCESS) {\r
449                         return;\r
450                 }\r
451                 \r
452                 // Create dense SUNMATRIX for use in linear solver\r
453                 // indexed by columns stacked on top of each other\r
454                 try {\r
455                     A = new double[NEQ][NEQ];\r
456                 }\r
457                 catch (OutOfMemoryError e) {\r
458                     MipavUtil.displayError("Out of memory error trying to create A");\r
459                     return;\r
460                 }\r
461                 \r
462                 // Create dense SUNLinearSolver object for use by CVode\r
463                 LS = SUNDenseLinearSolver(y, A);\r
464                 if (LS == null) {\r
465                         return;\r
466                 }\r
467                 \r
468                 // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
469                 //flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
470         }\r
471         \r
472         private void N_VNew_Serial(NVector y, int length) {\r
473             if ((y != null) && (length > 0)) {\r
474                 y.data = new double[length];\r
475                 y.own_data = true;\r
476             }\r
477         }\r
478         \r
479         /**\r
480          * f routine.  Compte function f(t,y).\r
481          * @param t\r
482          * @param yv\r
483          * @param ydotv\r
484          * @return\r
485          */\r
486         private int f(double t, NVector yv, NVector ydotv) {\r
487                 double y[] = yv.getData();\r
488                 double ydot[] = ydotv.getData();\r
489                 if (problem == 1) {\r
490                     ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
491                     ydot[2] = 3.0E7*y[1]*y[1];\r
492                     ydot[1] = -ydot[0] - ydot[2];\r
493                 }\r
494                 return 0;\r
495         }\r
496         \r
497         /**\r
498          * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
499          * @param t\r
500          * @param yv\r
501          * @param gout\r
502          * @return\r
503          */\r
504         private int g(double t, NVector yv, double gout[]) {\r
505                 double y[] = yv.getData();\r
506                 if (problem == 1) {\r
507                     gout[0] = y[0] - 0.0001;\r
508                     gout[1] = y[2] - 0.01;\r
509                 }\r
510                 return 0;\r
511         }\r
512         \r
513         // Types: struct CVodeMemRec, CVodeMem\r
514         // -----------------------------------------------------------------\r
515         // The type CVodeMem is type pointer to struct CVodeMemRec.\r
516         // This structure contains fields to keep track of problem state.\r
517         \r
518         private class NVector {\r
519                 double data[];\r
520                 boolean own_data;\r
521                 \r
522                 void setData(double data[]) {\r
523                         this.data = data;\r
524                 }\r
525                 double[] getData() {\r
526                     return data;        \r
527                 }\r
528         }\r
529         \r
530           \r
531         private class CVodeMemRec {\r
532             \r
533           double cv_uround;         /* machine unit roundoff                         */   \r
534 \r
535           /*-------------------------- \r
536             Problem Specification Data \r
537             --------------------------*/\r
538 \r
539           //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
540           int cv_f;\r
541           //void *cv_user_data;         /* user pointer passed to f                      */\r
542 \r
543           int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
544           int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
545 \r
546           int cv_itol;                /* itol = CV_SS, CV_SV, or CV_WF, or CV_NN       */\r
547           double cv_reltol;         /* relative tolerance                            */\r
548           double cv_Sabstol;        /* scalar absolute tolerance                     */\r
549           NVector cv_Vabstol = new NVector();        /* vector absolute tolerance                     */\r
550           boolean cv_user_efun;   /* SUNTRUE if user sets efun                     */\r
551           //CVEwtFn cv_efun;            /* function to set ewt                           */\r
552           //void *cv_e_data;            /* user pointer passed to efun                   */\r
553 \r
554           /*-----------------------\r
555             Quadrature Related Data \r
556             -----------------------*/\r
557 \r
558           boolean cv_quadr;       /* SUNTRUE if integrating quadratures            */\r
559 \r
560           //CVQuadRhsFn cv_fQ;          /* q' = fQ(t, y(t))                              */\r
561 \r
562           boolean cv_errconQ;     /* SUNTRUE if quadrs. are included in error test */\r
563 \r
564           int cv_itolQ;               /* itolQ = CV_SS or CV_SV                        */\r
565           double cv_reltolQ;        /* relative tolerance for quadratures            */\r
566           double cv_SabstolQ;       /* scalar absolute tolerance for quadratures     */\r
567           NVector cv_VabstolQ = new NVector();       /* vector absolute tolerance for quadratures     */\r
568 \r
569           /*------------------------\r
570             Sensitivity Related Data \r
571             ------------------------*/\r
572 \r
573           boolean cv_sensi;       /* SUNTRUE if computing sensitivities           */\r
574 \r
575           int cv_Ns;                  /* Number of sensitivities                      */\r
576 \r
577           int cv_ism;                 /* ism = SIMULTANEOUS or STAGGERED              */\r
578 \r
579           //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
580           //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
581           //void *cv_fS_data;           /* data pointer passed to fS                    */\r
582           boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
583           int cv_ifS;                 /* ifS = ALLSENS or ONESENS                     */\r
584 \r
585           double cv_p[];             /* parameters in f(t,y,p)                       */\r
586           double cv_pbar[];          /* scale factors for parameters                 */\r
587           int cv_plist[];              /* list of sensitivities                        */\r
588           int cv_DQtype;              /* central/forward finite differences           */\r
589           double cv_DQrhomax;       /* cut-off value for separate/simultaneous FD   */\r
590 \r
591           boolean cv_errconS;     /* SUNTRUE if yS are considered in err. control */\r
592 \r
593           int cv_itolS;\r
594           double cv_reltolS;        /* relative tolerance for sensitivities         */\r
595           double cv_SabstolS[];      /* scalar absolute tolerances for sensi.        */\r
596           NVector cv_VabstolS = new NVector();      /* vector absolute tolerances for sensi.        */\r
597 \r
598           /*-----------------------------------\r
599             Quadrature Sensitivity Related Data \r
600             -----------------------------------*/\r
601 \r
602           boolean cv_quadr_sensi; /* SUNTRUE if computing sensitivties of quadrs. */\r
603 \r
604           //CVQuadSensRhsFn cv_fQS;     /* fQS = (dfQ/dy)*yS + (dfQ/dp)                 */\r
605           //void *cv_fQS_data;          /* data pointer passed to fQS                   */\r
606           boolean cv_fQSDQ;       /* SUNTRUE if using internal DQ functions       */\r
607 \r
608           boolean cv_errconQS;    /* SUNTRUE if yQS are considered in err. con.   */\r
609 \r
610           int cv_itolQS;\r
611           double cv_reltolQS;       /* relative tolerance for yQS                   */\r
612           double cv_SabstolQS[];     /* scalar absolute tolerances for yQS           */\r
613           NVector cv_VabstolQS = new NVector();     /* vector absolute tolerances for yQS           */\r
614 \r
615           /*-----------------------\r
616             Nordsieck History Array \r
617             -----------------------*/\r
618 \r
619           NVector cv_zn[] = new NVector[L_MAX];   /* Nordsieck array, of size N x (q+1).\r
620                                          zn[j] is a vector of length N (j=0,...,q)\r
621                                          zn[j] = [1/factorial(j)] * h^j * \r
622                                          (jth derivative of the interpolating poly.)  */\r
623 \r
624           /*-------------------\r
625             Vectors of length N \r
626             -------------------*/\r
627 \r
628           NVector cv_ewt = new NVector();            /* error weight vector                          */\r
629           NVector cv_y = new NVector();              // y is used as temporary storage by the solver.\r
630                                       /*   The memory is provided by the user to CVode \r
631                                          where the vector is named yout.              */\r
632           NVector cv_acor = new NVector();           // In the context of the solution of the\r
633                                        /*  nonlinear equation, acor = y_n(m) - y_n(0).\r
634                                          On return, this vector is scaled to give\r
635                                          the estimated local error in y.              */\r
636           NVector cv_tempv = new NVector();          /* temporary storage vector                     */\r
637           NVector cv_ftemp = new NVector();          /* temporary storage vector                     */\r
638           \r
639           /*--------------------------\r
640             Quadrature Related Vectors \r
641             --------------------------*/\r
642 \r
643           NVector cv_znQ[] = new NVector[L_MAX];     /* Nordsieck arrays for quadratures             */\r
644           NVector cv_ewtQ = new NVector();           /* error weight vector for quadratures          */\r
645           NVector cv_yQ = new NVector();             /* Unlike y, yQ is not allocated by the user    */\r
646           NVector cv_acorQ = new NVector();          /* acorQ = yQ_n(m) - yQ_n(0)                    */\r
647           NVector cv_tempvQ = new NVector();         /* temporary storage vector (~ tempv)           */\r
648 \r
649           /*---------------------------\r
650             Sensitivity Related Vectors \r
651             ---------------------------*/\r
652 \r
653           NVector cv_znS[] = new NVector[L_MAX];    /* Nordsieck arrays for sensitivities           */\r
654           NVector cv_ewtS[];          /* error weight vectors for sensitivities       */\r
655           NVector cv_yS[];            /* yS=yS0 (allocated by the user)               */\r
656           NVector cv_acorS[];         /* acorS = yS_n(m) - yS_n(0)                    */\r
657           NVector cv_tempvS[];        /* temporary storage vector (~ tempv)           */\r
658           NVector cv_ftempS[];        /* temporary storage vector (~ ftemp)           */\r
659 \r
660           boolean cv_stgr1alloc;  /* Did we allocate ncfS1, ncfnS1, and nniS1?    */\r
661 \r
662           /*--------------------------------------\r
663             Quadrature Sensitivity Related Vectors \r
664             --------------------------------------*/\r
665 \r
666           NVector cv_znQS[] = new NVector[L_MAX];   /* Nordsieck arrays for quadr. sensitivities    */\r
667           NVector cv_ewtQS[];         /* error weight vectors for sensitivities       */\r
668           NVector cv_yQS[];           /* Unlike yS, yQS is not allocated by the user  */\r
669           NVector cv_acorQS[];        /* acorQS = yQS_n(m) - yQS_n(0)                 */\r
670           NVector cv_tempvQS[];       /* temporary storage vector (~ tempv)           */\r
671           NVector cv_ftempQ = new NVector();         /* temporary storage vector (~ ftemp)           */\r
672           \r
673           /*-----------------\r
674             Tstop information\r
675             -----------------*/\r
676 \r
677           boolean cv_tstopset;\r
678           double cv_tstop;\r
679 \r
680           /*---------\r
681             Step Data \r
682             ---------*/\r
683 \r
684           int cv_q;                    /* current order                               */\r
685           int cv_qprime;               /* order to be used on the next step\r
686                                         * qprime = q-1, q, or q+1                     */\r
687           int cv_next_q;               /* order to be used on the next step           */\r
688           int cv_qwait;                /* number of internal steps to wait before\r
689                                         * considering a change in q                   */\r
690           int cv_L;                    /* L = q + 1                                   */\r
691 \r
692           double cv_hin;\r
693           double cv_h;               /* current step size                           */\r
694           double cv_hprime;          /* step size to be used on the next step       */ \r
695           double cv_next_h;          /* step size to be used on the next step       */ \r
696           double cv_eta;             /* eta = hprime / h                            */\r
697           double cv_hscale;          /* value of h used in zn                       */\r
698           double cv_tn;              /* current internal value of t                 */\r
699           double cv_tretlast;        /* last value of t returned                    */\r
700 \r
701           double cv_tau[] = new double[L_MAX+1];    /* array of previous q+1 successful step\r
702                                         * sizes indexed from 1 to q+1                 */\r
703           double cv_tq[] =new double[NUM_TESTS+1]; /* array of test quantities indexed from\r
704                                         * 1 to NUM_TESTS(=5)                          */\r
705           double cv_l[] = new double[L_MAX];        /* coefficients of l(x) (degree q poly)        */\r
706 \r
707           double cv_rl1;             /* the scalar 1/l[1]                           */\r
708           double cv_gamma;           /* gamma = h * rl1                             */\r
709           double cv_gammap;          /* gamma at the last setup call                */\r
710           double cv_gamrat;          /* gamma / gammap                              */\r
711 \r
712           double cv_crate;           /* est. corrector conv. rate in Nls            */\r
713           double cv_crateS;          /* est. corrector conv. rate in NlsStgr        */\r
714           double cv_acnrm;           /* | acor |                                    */\r
715           double cv_acnrmQ;          /* | acorQ |                                   */\r
716           double cv_acnrmS;          /* | acorS |                                   */\r
717           double cv_acnrmQS;         /* | acorQS |                                  */\r
718           double cv_nlscoef;         /* coeficient in nonlinear convergence test    */\r
719           int  cv_mnewt;               /* Newton iteration counter                    */\r
720           int  cv_ncfS1[];              /* Array of Ns local counters for conv.  \r
721                                         * failures (used in CVStep for STAGGERED1)    */\r
722 \r
723           /*------\r
724             Limits \r
725             ------*/\r
726 \r
727           int cv_qmax;             /* q <= qmax                                       */\r
728           long cv_mxstep;      /* maximum number of internal steps for one \r
729                                       user call                                       */\r
730           int cv_maxcor;           /* maximum number of corrector iterations for \r
731                                       the solution of the nonlinear equation          */\r
732           int cv_maxcorS;\r
733           int cv_mxhnil;           /* max. number of warning messages issued to the\r
734                                       user that t + h == t for the next internal step */\r
735           int cv_maxnef;           /* maximum number of error test failures           */\r
736           int cv_maxncf;           /* maximum number of nonlinear conv. failures      */\r
737           \r
738           double cv_hmin;        /* |h| >= hmin                                     */\r
739           double cv_hmax_inv;    /* |h| <= 1/hmax_inv                               */\r
740           double cv_etamax;      /* eta <= etamax                                   */\r
741 \r
742           /*----------\r
743             Counters \r
744             ----------*/\r
745 \r
746           long cv_nst;         /* number of internal steps taken                  */\r
747 \r
748           long cv_nfe;         /* number of f calls                               */\r
749           long cv_nfQe;        /* number of fQ calls                              */\r
750           long cv_nfSe;        /* number of fS calls                              */\r
751           long cv_nfeS;        /* number of f calls from sensi DQ                 */\r
752           long cv_nfQSe;       /* number of fQS calls                             */\r
753           long cv_nfQeS;       /* number of fQ calls from sensi DQ                */\r
754 \r
755 \r
756           long cv_ncfn;        /* number of corrector convergence failures        */\r
757           long cv_ncfnS;       /* number of total sensi. corr. conv. failures     */\r
758           long cv_ncfnS1[];     /* number of sensi. corrector conv. failures       */\r
759 \r
760           long cv_nni;         /* number of nonlinear iterations performed        */\r
761           long cv_nniS;        /* number of total sensi. nonlinear iterations     */\r
762           long cv_nniS1[];      /* number of sensi. nonlinear iterations           */\r
763 \r
764           long cv_netf;        /* number of error test failures                   */\r
765           long cv_netfQ;       /* number of quadr. error test failures            */\r
766           long cv_netfS;       /* number of sensi. error test failures            */\r
767           long cv_netfQS;      /* number of quadr. sensi. error test failures     */\r
768 \r
769           long cv_nsetups;     /* number of setup calls                           */\r
770           long cv_nsetupsS;    /* number of setup calls due to sensitivities      */\r
771 \r
772           int cv_nhnil;            /* number of messages issued to the user that\r
773                                       t + h == t for the next iternal step            */\r
774 \r
775           /*-----------------------------\r
776             Space requirements for CVODES \r
777             -----------------------------*/\r
778 \r
779           long cv_lrw1;        /* no. of realtype words in 1 N_Vector y           */ \r
780           long cv_liw1;        /* no. of integer words in 1 N_Vector y            */ \r
781           long cv_lrw1Q;       /* no. of realtype words in 1 N_Vector yQ          */ \r
782           long cv_liw1Q;       /* no. of integer words in 1 N_Vector yQ           */ \r
783           long cv_lrw;             /* no. of realtype words in CVODES work vectors    */\r
784           long cv_liw;             /* no. of integer words in CVODES work vectors     */\r
785 \r
786           /*----------------\r
787             Step size ratios\r
788             ----------------*/\r
789 \r
790           double cv_etaqm1;      /* ratio of new to old h for order q-1             */\r
791           double cv_etaq;        /* ratio of new to old h for order q               */\r
792           double cv_etaqp1;      /* ratio of new to old h for order q+1             */\r
793 \r
794           /*------------------\r
795             Linear Solver Data \r
796             ------------------*/\r
797 \r
798           /* Linear Solver functions to be called */\r
799 \r
800           //int (*cv_linit)(struct CVodeMemRec *cv_mem);\r
801 \r
802           //int (*cv_lsetup)(struct CVodeMemRec *cv_mem, int convfail, \r
803                            //N_Vector ypred, N_Vector fpred, booleantype *jcurPtr, \r
804                            //N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3); \r
805 \r
806           //int (*cv_lsolve)(struct CVodeMemRec *cv_mem, N_Vector b, N_Vector weight,\r
807                            //N_Vector ycur, N_Vector fcur);\r
808 \r
809           //int (*cv_lfree)(struct CVodeMemRec *cv_mem);\r
810 \r
811           /* Linear Solver specific memory */\r
812 \r
813           //void *cv_lmem;           \r
814 \r
815           /* Flag to request a call to the setup routine */\r
816 \r
817           boolean cv_forceSetup;\r
818 \r
819           /*------------\r
820             Saved Values\r
821             ------------*/\r
822 \r
823           int cv_qu;                   /* last successful q value used                */\r
824           long cv_nstlp;           /* step number of last setup call              */\r
825           double cv_h0u;             /* actual initial stepsize                     */\r
826           double cv_hu;              /* last successful h value used                */\r
827           double cv_saved_tq5;       /* saved value of tq[5]                        */\r
828           boolean cv_jcur;         /* is Jacobian info for linear solver current? */\r
829           int cv_convfail;             /* flag storing previous solver failure mode   */\r
830           double cv_tolsf;           /* tolerance scale factor                      */\r
831           int cv_qmax_alloc;           /* qmax used when allocating mem               */\r
832           int cv_qmax_allocQ;          /* qmax used when allocating quad. mem         */\r
833           int cv_qmax_allocS;          /* qmax used when allocating sensi. mem        */\r
834           int cv_qmax_allocQS;         /* qmax used when allocating quad. sensi. mem  */\r
835           int cv_indx_acor;            /* index of zn vector in which acor is saved   */\r
836 \r
837           /*--------------------------------------------------------------------\r
838             Flags turned ON by CVodeInit, CVodeSensMalloc, and CVodeQuadMalloc \r
839             and read by CVodeReInit, CVodeSensReInit, and CVodeQuadReInit\r
840             --------------------------------------------------------------------*/\r
841 \r
842           boolean cv_VabstolMallocDone;\r
843           boolean cv_MallocDone;\r
844 \r
845           boolean cv_VabstolQMallocDone;\r
846           boolean cv_QuadMallocDone;\r
847 \r
848           boolean cv_VabstolSMallocDone;\r
849           boolean cv_SabstolSMallocDone;\r
850           boolean cv_SensMallocDone;\r
851 \r
852           boolean cv_VabstolQSMallocDone;\r
853           boolean cv_SabstolQSMallocDone;\r
854           boolean cv_QuadSensMallocDone;\r
855 \r
856           /*-------------------------------------------\r
857             Error handler function and error ouput file \r
858             -------------------------------------------*/\r
859 \r
860           //CVErrHandlerFn cv_ehfun;    /* Error messages are handled by ehfun          */\r
861           CVodeMemRec cv_eh_data;           /* dats pointer passed to ehfun                 */\r
862           RandomAccessFile cv_errfp;             /* CVODES error messages are sent to errfp      */    \r
863 \r
864           /*-------------------------\r
865             Stability Limit Detection\r
866             -------------------------*/\r
867 \r
868           boolean cv_sldeton;     /* Is Stability Limit Detection on?             */\r
869           double cv_ssdat[][] = new double[6][4];    /* scaled data array for STALD                  */\r
870           int cv_nscon;               /* counter for STALD method                     */\r
871           long cv_nor;            /* counter for number of order reductions       */\r
872 \r
873           /*----------------\r
874             Rootfinding Data\r
875             ----------------*/\r
876 \r
877           //CVRootFn cv_gfun;        /* Function g for roots sought                     */\r
878           int cv_gfun;\r
879           int cv_nrtfn;            /* number of components of g                       */\r
880           int cv_iroots[];          /* array for root information                      */\r
881           int cv_rootdir[];         /* array specifying direction of zero-crossing     */\r
882           double cv_tlo;         /* nearest endpoint of interval in root search     */\r
883           double cv_thi;         /* farthest endpoint of interval in root search    */\r
884           double cv_trout;       /* t value returned by rootfinding routine         */\r
885           double cv_glo[];        /* saved array of g values at t = tlo              */\r
886           double cv_ghi[];        /* saved array of g values at t = thi              */\r
887           double cv_grout[];      /* array of g values at t = trout                  */\r
888           double cv_toutc;       /* copy of tout (if NORMAL mode)                   */\r
889           double cv_ttol;        /* tolerance on root location trout                */\r
890           int cv_taskc;            /* copy of parameter itask                         */\r
891           int cv_irfnd;            /* flag showing whether last step had a root       */\r
892           long cv_nge;         /* counter for g evaluations                       */\r
893           boolean cv_gactive[]; /* array with active/inactive event functions      */\r
894           int cv_mxgnull;          /* number of warning messages about possible g==0  */\r
895 \r
896           /*------------------------\r
897             Adjoint sensitivity data\r
898             ------------------------*/\r
899 \r
900           boolean cv_adj;             /* SUNTRUE if performing ASA                */\r
901 \r
902           CVodeMemRec cv_adj_mem; /* Pointer to adjoint memory structure      */\r
903 \r
904           boolean cv_adjMallocDone;\r
905 \r
906         } ;\r
907 \r
908 \r
909         \r
910         // Call CVodeCreate to create the solver memory and specify the\r
911         // Backward Differentiation Formula and the use of a Newton \r
912         // iteration\r
913         \r
914          // CVodeCreate creates an internal memory block for a problem to\r
915          // be solved by CVODES.\r
916          //\r
917          // lmm  is the type of linear multistep method to be used.\r
918          //      The legal values are CV_ADAMS and CV_BDF (see previous\r
919          //      description).\r
920          //\r
921          // iter  is the type of iteration used to solve the nonlinear\r
922          //       system that arises during each internal time step.\r
923          //       The legal values are CV_FUNCTIONAL and CV_NEWTON.\r
924          //\r
925          //If successful, CVodeCreate returns a pointer to initialized\r
926          // problem memory. This pointer should be passed to CVodeInit.\r
927          // If an initialization error occurs, CVodeCreate prints an error\r
928          // message to standard err and returns NULL.\r
929 \r
930      private CVodeMemRec CVodeCreate(int lmm, int iter) {\r
931          int maxord;\r
932          CVodeMemRec cv_mem;\r
933 \r
934           /* Test inputs */\r
935 \r
936           if ((lmm != CV_ADAMS) && (lmm != CV_BDF)) {\r
937             cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_LMM);\r
938             return null;\r
939           }\r
940           \r
941           if ((iter != CV_FUNCTIONAL) && (iter != CV_NEWTON)) {\r
942             cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_ITER);\r
943             return null;\r
944           }\r
945 \r
946           cv_mem = null;\r
947           cv_mem = new CVodeMemRec();\r
948 \r
949           /* Zero out cv_mem */\r
950           //memset(cv_mem, 0, sizeof(struct CVodeMemRec));\r
951 \r
952           maxord = (lmm == CV_ADAMS) ? ADAMS_Q_MAX : BDF_Q_MAX;\r
953 \r
954           /* copy input parameters into cv_mem */\r
955 \r
956           cv_mem.cv_lmm  = lmm;\r
957           cv_mem.cv_iter = iter;\r
958 \r
959           /* Set uround */\r
960 \r
961           cv_mem.cv_uround = UNIT_ROUNDOFF;\r
962 \r
963           /* Set default values for integrator optional inputs */\r
964 \r
965           //cv_mem.cv_f          = null;\r
966           cv_mem.cv_f = -1;\r
967           //cv_mem.cv_user_data  = null;\r
968           cv_mem.cv_itol       = CV_NN;\r
969           cv_mem.cv_user_efun  = false;\r
970           //cv_mem.cv_efun       = null;\r
971           //cv_mem.cv_e_data     = null;\r
972           //cv_mem.cv_ehfun      = cvErrHandler;\r
973           cv_mem.cv_eh_data    = cv_mem;\r
974           //cv_mem.cv_errfp      = stderr;\r
975           cv_mem.cv_qmax       = maxord;\r
976           cv_mem.cv_mxstep     = MXSTEP_DEFAULT;\r
977           cv_mem.cv_mxhnil     = MXHNIL_DEFAULT;\r
978           cv_mem.cv_sldeton    = false;\r
979           cv_mem.cv_hin        = ZERO;\r
980           cv_mem.cv_hmin       = HMIN_DEFAULT;\r
981           cv_mem.cv_hmax_inv   = HMAX_INV_DEFAULT;\r
982           cv_mem.cv_tstopset   = false;\r
983           cv_mem.cv_maxcor     = NLS_MAXCOR;\r
984           cv_mem.cv_maxnef     = MXNEF;\r
985           cv_mem.cv_maxncf     = MXNCF;\r
986           cv_mem.cv_nlscoef    = CORTES;\r
987 \r
988           /* Initialize root finding variables */\r
989 \r
990           cv_mem.cv_glo        = null;\r
991           cv_mem.cv_ghi        = null;\r
992           cv_mem.cv_grout      = null;\r
993           cv_mem.cv_iroots     = null;\r
994           cv_mem.cv_rootdir    = null;\r
995           //cv_mem.cv_gfun       = null;\r
996           cv_mem.cv_gfun = -1;\r
997           cv_mem.cv_nrtfn      = 0;  \r
998           cv_mem.cv_gactive    = null;\r
999           cv_mem.cv_mxgnull    = 1;\r
1000 \r
1001           /* Set default values for quad. optional inputs */\r
1002 \r
1003           cv_mem.cv_quadr      = false;\r
1004           //cv_mem.cv_fQ         = null;\r
1005           cv_mem.cv_errconQ    = false;\r
1006           cv_mem.cv_itolQ      = CV_NN;\r
1007 \r
1008           /* Set default values for sensi. optional inputs */\r
1009 \r
1010           cv_mem.cv_sensi      = false;\r
1011           //cv_mem.cv_fS_data    = mull;\r
1012           //cv_mem.cv_fS         = cvSensRhsInternalDQ;\r
1013           //cv_mem.cv_fS1        = cvSensRhs1InternalDQ;\r
1014           cv_mem.cv_fSDQ       = true;\r
1015           cv_mem.cv_ifS        = CV_ONESENS;\r
1016           cv_mem.cv_DQtype     = CV_CENTERED;\r
1017           cv_mem.cv_DQrhomax   = ZERO;\r
1018           cv_mem.cv_p          = null;\r
1019           cv_mem.cv_pbar       = null;\r
1020           cv_mem.cv_plist      = null;\r
1021           cv_mem.cv_errconS    = false;\r
1022           cv_mem.cv_maxcorS    = NLS_MAXCOR;\r
1023           cv_mem.cv_ncfS1      = null;\r
1024           cv_mem.cv_ncfnS1     = null;\r
1025           cv_mem.cv_nniS1      = null;\r
1026           cv_mem.cv_itolS      = CV_NN;\r
1027 \r
1028           /* Set default values for quad. sensi. optional inputs */\r
1029 \r
1030           cv_mem.cv_quadr_sensi = false;\r
1031           //cv_mem.cv_fQS         = null;\r
1032           //cv_mem.cv_fQS_data    = null;\r
1033           cv_mem.cv_fQSDQ       = true;\r
1034           cv_mem.cv_errconQS    = false;\r
1035           cv_mem.cv_itolQS      = CV_NN;\r
1036 \r
1037           /* Set default for ASA */\r
1038 \r
1039           cv_mem.cv_adj         = false;\r
1040           cv_mem.cv_adj_mem     = null;\r
1041 \r
1042           /* Set the saved values for qmax_alloc */\r
1043 \r
1044           cv_mem.cv_qmax_alloc  = maxord;\r
1045           cv_mem.cv_qmax_allocQ = maxord;\r
1046           cv_mem.cv_qmax_allocS = maxord;\r
1047 \r
1048           /* Initialize lrw and liw */\r
1049 \r
1050           cv_mem.cv_lrw = 65 + 2*L_MAX + NUM_TESTS;\r
1051           cv_mem.cv_liw = 52;\r
1052 \r
1053           /* No mallocs have been done yet */\r
1054 \r
1055           cv_mem.cv_VabstolMallocDone   = false;\r
1056           cv_mem.cv_MallocDone          = false;\r
1057 \r
1058           cv_mem.cv_VabstolQMallocDone  = false;\r
1059           cv_mem.cv_QuadMallocDone      = false;\r
1060 \r
1061           cv_mem.cv_VabstolSMallocDone  = false;\r
1062           cv_mem.cv_SabstolSMallocDone  = false;\r
1063           cv_mem.cv_SensMallocDone      = false;\r
1064 \r
1065           cv_mem.cv_VabstolQSMallocDone = false;\r
1066           cv_mem.cv_SabstolQSMallocDone = false;\r
1067           cv_mem.cv_QuadSensMallocDone  = false;\r
1068 \r
1069           cv_mem.cv_adjMallocDone       = false;\r
1070 \r
1071           /* Return pointer to CVODES memory block */\r
1072 \r
1073           return cv_mem;\r
1074 \r
1075      }\r
1076      \r
1077      /*\r
1078       * cvProcessError is a high level error handling function.\r
1079       * - If cv_mem==NULL it prints the error message to stderr.\r
1080       * - Otherwise, it sets up and calls the error handling function \r
1081       *   pointed to by cv_ehfun.\r
1082       */\r
1083 \r
1084      void cvProcessError(CVodeMemRec cv_mem, \r
1085                          int error_code,  String module, String fname, \r
1086                          String msgfmt)\r
1087      {\r
1088          MipavUtil.displayError("In module " + module + " in function " + fname + " msgfmt");\r
1089        //va_list ap;\r
1090        //char msg[256];\r
1091 \r
1092        /* Initialize the argument pointer variable \r
1093           (msgfmt is the last required argument to cvProcessError) */\r
1094 \r
1095        //va_start(ap, msgfmt);\r
1096 \r
1097        /* Compose the message */\r
1098 \r
1099        //vsprintf(msg, msgfmt, ap);\r
1100 \r
1101        //if (cv_mem == NULL) {    /* We write to stderr */\r
1102      //#ifndef NO_FPRINTF_OUTPUT\r
1103          //fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);\r
1104          //fprintf(stderr, "%s\n\n", msg);\r
1105      //#endif\r
1106 \r
1107        //} else {                 /* We can call ehfun */\r
1108          //cv_mem->cv_ehfun(error_code, module, fname, msg, cv_mem->cv_eh_data);\r
1109        //}\r
1110 \r
1111        /* Finalize argument processing */\r
1112        //va_end(ap);\r
1113 \r
1114        return;\r
1115      }\r
1116      \r
1117      /*-----------------------------------------------------------------*/\r
1118 \r
1119      /*\r
1120       * CVodeInit\r
1121       * \r
1122       * CVodeInit allocates and initializes memory for a problem. All \r
1123       * problem inputs are checked for errors. If any error occurs during \r
1124       * initialization, it is reported to the file whose file pointer is \r
1125       * errfp and an error flag is returned. Otherwise, it returns CV_SUCCESS\r
1126       */\r
1127 \r
1128      int CVodeInit(CVodeMemRec cv_mem, int f, double t0, NVector y0)\r
1129      {\r
1130        boolean nvectorOK, allocOK;\r
1131        long lrw1[] = new long[1];\r
1132        long liw1[] = new long[1];\r
1133        int i,k;\r
1134 \r
1135        /* Check cvode_mem */\r
1136 \r
1137        if (cv_mem==null) {\r
1138          cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeInit",\r
1139                         MSGCV_NO_MEM);\r
1140          return(CV_MEM_NULL);\r
1141        }\r
1142 \r
1143        /* Check for legal input parameters */\r
1144 \r
1145        if (y0==null) {\r
1146          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1147                         MSGCV_NULL_Y0);\r
1148          return(CV_ILL_INPUT);\r
1149        }\r
1150 \r
1151        if (f < 0) {\r
1152          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1153                         MSGCV_NULL_F);\r
1154          return(CV_ILL_INPUT);\r
1155        }\r
1156 \r
1157        /* Test if all required vector operations are implemented */\r
1158 \r
1159        /*nvectorOK = cvCheckNvector(y0);\r
1160        if(!nvectorOK) {\r
1161          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1162                         MSGCV_BAD_NVECTOR);\r
1163          return(CV_ILL_INPUT);\r
1164        }*/\r
1165 \r
1166        /* Set space requirements for one N_Vector */\r
1167 \r
1168        if (y0.getData() != null) {\r
1169          N_VSpace_Serial(y0, lrw1, liw1);\r
1170        } else {\r
1171          lrw1[0] = 0;\r
1172          liw1[0] = 0;\r
1173        }\r
1174        cv_mem.cv_lrw1 = lrw1[0];\r
1175        cv_mem.cv_liw1 = liw1[0];\r
1176 \r
1177        /* Allocate the vectors (using y0 as a template) */\r
1178 \r
1179       allocOK = cvAllocVectors(cv_mem, y0);\r
1180        if (!allocOK) {\r
1181          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeInit",\r
1182                         MSGCV_MEM_FAIL);\r
1183          return(CV_MEM_FAIL);\r
1184        }\r
1185 \r
1186        /* All error checking is complete at this point */\r
1187 \r
1188        /* Copy the input parameters into CVODES state */\r
1189 \r
1190        //cv_mem.cv_f  = f;\r
1191        cv_mem.cv_tn = t0;\r
1192 \r
1193        /* Set step parameters */\r
1194 \r
1195        cv_mem.cv_q      = 1;\r
1196        cv_mem.cv_L      = 2;\r
1197        cv_mem.cv_qwait  = cv_mem.cv_L;\r
1198        cv_mem.cv_etamax = ETAMX1;\r
1199 \r
1200        cv_mem.cv_qu     = 0;\r
1201        cv_mem.cv_hu     = ZERO;\r
1202        cv_mem.cv_tolsf  = ONE;\r
1203 \r
1204        /* Set the linear solver addresses to NULL.\r
1205           (We check != NULL later, in CVode, if using CV_NEWTON.) */\r
1206 \r
1207        //cv_mem.cv_linit  = null;\r
1208        //cv_mem.cv_lsetup = null;\r
1209        //cv_mem.cv_lsolve = null;\r
1210        //cv_mem.cv_lfree  = null;\r
1211        //cv_mem.cv_lmem   = null;\r
1212 \r
1213        /* Set forceSetup to SUNFALSE */\r
1214 \r
1215        cv_mem.cv_forceSetup = false;\r
1216 \r
1217        /* Initialize zn[0] in the history array */\r
1218 \r
1219        N_VScale(ONE, y0, cv_mem.cv_zn[0]);\r
1220 \r
1221        /* Initialize all the counters */\r
1222 \r
1223        cv_mem.cv_nst     = 0;\r
1224        cv_mem.cv_nfe     = 0;\r
1225        cv_mem.cv_ncfn    = 0;\r
1226        cv_mem.cv_netf    = 0;\r
1227        cv_mem.cv_nni     = 0;\r
1228        cv_mem.cv_nsetups = 0;\r
1229        cv_mem.cv_nhnil   = 0;\r
1230        cv_mem.cv_nstlp   = 0;\r
1231        cv_mem.cv_nscon   = 0;\r
1232        cv_mem.cv_nge     = 0;\r
1233 \r
1234        cv_mem.cv_irfnd   = 0;\r
1235 \r
1236        /* Initialize other integrator optional outputs */\r
1237 \r
1238        cv_mem.cv_h0u      = ZERO;\r
1239        cv_mem.cv_next_h   = ZERO;\r
1240        cv_mem.cv_next_q   = 0;\r
1241 \r
1242        /* Initialize Stablilty Limit Detection data */\r
1243        /* NOTE: We do this even if stab lim det was not\r
1244           turned on yet. This way, the user can turn it\r
1245           on at any time */\r
1246 \r
1247        cv_mem.cv_nor = 0;\r
1248        for (i = 1; i <= 5; i++)\r
1249          for (k = 1; k <= 3; k++) \r
1250            cv_mem.cv_ssdat[i-1][k-1] = ZERO;\r
1251 \r
1252        /* Problem has been successfully initialized */\r
1253 \r
1254        cv_mem.cv_MallocDone = true;\r
1255 \r
1256        return(CV_SUCCESS);\r
1257      }\r
1258      \r
1259      void N_VSpace_Serial(NVector v, long lrw[], long liw[])\r
1260      {\r
1261        lrw[0] = v.getData().length;\r
1262        liw[0] = 1;\r
1263 \r
1264        return;\r
1265      }\r
1266 \r
1267      \r
1268      /*\r
1269       * cvCheckNvector\r
1270       * This routine checks if all required vector operations are present.\r
1271       * If any of them is missing it returns SUNFALSE.\r
1272       */\r
1273 \r
1274      /*static booleantype cvCheckNvector(N_Vector tmpl)\r
1275      {\r
1276        if((tmpl->ops->nvclone     == NULL) || // clone\r
1277           (tmpl->ops->nvdestroy   == NULL) || // set null\r
1278           (tmpl->ops->nvlinearsum == NULL) || // (double a, double x[], double b, double y[], double z[])\r
1279                                               // z[i] = a*x[i] + b*y[i];\r
1280           (tmpl->ops->nvconst     == NULL) || // set all to constant\r
1281           (tmpl->ops->nvprod      == NULL) || // (double x[], double y[], double z[])\r
1282                                               // z[i] = x[i] * y[i]\r
1283           (tmpl->ops->nvdiv       == NULL) || // double x[], double y[], double z[])\r
1284                                               // z[i] = x[i]/y[i]\r
1285           (tmpl->ops->nvscale     == NULL) || // (double c, double x[], double z[])\r
1286                                               // z[i] = c * x[i]\r
1287           (tmpl->ops->nvabs       == NULL) || // (double x[], double z[])\r
1288                                               // z[i] = abs(x[i])\r
1289           (tmpl->ops->nvinv       == NULL) || // (double z[], double z[])\r
1290                                               // z[i] = 1.0/x[i]\r
1291           (tmpl->ops->nvaddconst  == NULL) || // (double x[], double b, double z[])\r
1292                                               // z[i] = x[i] + b\r
1293           (tmpl->ops->nvmaxnorm   == NULL) || // (double x[])\r
1294                                               // max(abs(x[i]))\r
1295           (tmpl->ops->nvwrmsnorm  == NULL) || // (double x[], double w[])\r
1296                                               // prod = x[i] * w[i]\r
1297                                               // sum = sum over all i of prod[i]*prod[i]\r
1298                                               // answer = sqrt(sum/N)\r
1299           (tmpl->ops->nvmin       == NULL))   // (double x[])\r
1300                                               // min(x[i])\r
1301          return(SUNFALSE);\r
1302        else\r
1303          return(SUNTRUE);\r
1304      }*/\r
1305      \r
1306      /*\r
1307       * cvAllocVectors\r
1308       *\r
1309       * This routine allocates the CVODES vectors ewt, acor, tempv, ftemp, and\r
1310       * zn[0], ..., zn[maxord].\r
1311       * If all memory allocations are successful, cvAllocVectors returns SUNTRUE. \r
1312       * Otherwise all allocated memory is freed and cvAllocVectors returns SUNFALSE.\r
1313       * This routine also sets the optional outputs lrw and liw, which are\r
1314       * (respectively) the lengths of the real and integer work spaces\r
1315       * allocated here.\r
1316       */\r
1317 \r
1318      private boolean cvAllocVectors(CVodeMemRec cv_mem, NVector tmpl)\r
1319      {\r
1320        int i, j;\r
1321 \r
1322        /* Allocate ewt, acor, tempv, ftemp */\r
1323        \r
1324        cv_mem.cv_ewt = N_VClone(tmpl);\r
1325        if (cv_mem.cv_ewt == null) return(false);\r
1326 \r
1327        cv_mem.cv_acor = N_VClone(tmpl);\r
1328        if (cv_mem.cv_acor == null) {\r
1329          N_VDestroy(cv_mem.cv_ewt);\r
1330          return false;\r
1331        }\r
1332 \r
1333        cv_mem.cv_tempv = N_VClone(tmpl);\r
1334        if (cv_mem.cv_tempv == null) {\r
1335          N_VDestroy(cv_mem.cv_ewt);\r
1336          N_VDestroy(cv_mem.cv_acor);\r
1337          return false;\r
1338        }\r
1339 \r
1340        cv_mem.cv_ftemp = N_VClone(tmpl);\r
1341        if (cv_mem.cv_ftemp == null) {\r
1342          N_VDestroy(cv_mem.cv_tempv);\r
1343          N_VDestroy(cv_mem.cv_ewt);\r
1344          N_VDestroy(cv_mem.cv_acor);\r
1345          return false;\r
1346        }\r
1347 \r
1348        /* Allocate zn[0] ... zn[qmax] */\r
1349 \r
1350        for (j=0; j <= cv_mem.cv_qmax; j++) {\r
1351          cv_mem.cv_zn[j] = N_VClone(tmpl);\r
1352          if (cv_mem.cv_zn[j] == null) {\r
1353            N_VDestroy(cv_mem.cv_ewt);\r
1354            N_VDestroy(cv_mem.cv_acor);\r
1355            N_VDestroy(cv_mem.cv_tempv);\r
1356            N_VDestroy(cv_mem.cv_ftemp);\r
1357            for (i=0; i < j; i++) N_VDestroy(cv_mem.cv_zn[i]);\r
1358            return false;\r
1359          }\r
1360        }\r
1361 \r
1362        /* Update solver workspace lengths  */\r
1363        cv_mem.cv_lrw += (cv_mem.cv_qmax + 5)*cv_mem.cv_lrw1;\r
1364        cv_mem.cv_liw += (cv_mem.cv_qmax + 5)*cv_mem.cv_liw1;\r
1365 \r
1366        /* Store the value of qmax used here */\r
1367        cv_mem.cv_qmax_alloc = cv_mem.cv_qmax;\r
1368 \r
1369        return true;\r
1370      }\r
1371      \r
1372      private NVector N_VClone(NVector y) {\r
1373          NVector x = new NVector();\r
1374          x.data = y.data.clone();\r
1375          x.own_data = y.own_data;\r
1376          return x;\r
1377      }\r
1378 \r
1379      private void N_VDestroy(NVector x) {\r
1380          x.data = null;\r
1381          x = null;\r
1382      }\r
1383      \r
1384      private void N_VScale(double c, NVector x, NVector z) {\r
1385          int i;\r
1386          double xa[] = x.getData();\r
1387          int n = xa.length;\r
1388          double za[] = z.getData();\r
1389          for (i = 0; i < n; i++) {\r
1390                  za[i] = c * xa[i];\r
1391          }\r
1392      }\r
1393      \r
1394      int CVodeSVtolerances(CVodeMemRec cv_mem, double reltol, NVector abstol)\r
1395      {\r
1396 \r
1397        if (cv_mem==null) {\r
1398          cvProcessError(null, CV_MEM_NULL, "CVODES",\r
1399                         "CVodeSVtolerances", MSGCV_NO_MEM);\r
1400          return(CV_MEM_NULL);\r
1401        }\r
1402 \r
1403        if (cv_mem.cv_MallocDone == false) {\r
1404          cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES",\r
1405                         "CVodeSVtolerances", MSGCV_NO_MALLOC);\r
1406          return(CV_NO_MALLOC);\r
1407        }\r
1408 \r
1409        /* Check inputs */\r
1410 \r
1411        if (reltol < ZERO) {\r
1412          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
1413                         "CVodeSVtolerances", MSGCV_BAD_RELTOL);\r
1414          return(CV_ILL_INPUT);\r
1415        }\r
1416 \r
1417        if (N_VMin(abstol) < ZERO) {\r
1418          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
1419                         "CVodeSVtolerances", MSGCV_BAD_ABSTOL);\r
1420          return(CV_ILL_INPUT);\r
1421        }\r
1422 \r
1423        /* Copy tolerances into memory */\r
1424        \r
1425        if ( !(cv_mem.cv_VabstolMallocDone) ) {\r
1426          cv_mem.cv_Vabstol = N_VClone(cv_mem.cv_ewt);\r
1427          cv_mem.cv_lrw += cv_mem.cv_lrw1;\r
1428          cv_mem.cv_liw += cv_mem.cv_liw1;\r
1429          cv_mem.cv_VabstolMallocDone = true;\r
1430        }\r
1431 \r
1432        cv_mem.cv_reltol = reltol;\r
1433        N_VScale(ONE, abstol, cv_mem.cv_Vabstol);\r
1434 \r
1435        cv_mem.cv_itol = CV_SV;\r
1436 \r
1437        cv_mem.cv_user_efun = false;\r
1438        //cv_mem.cv_efun = cvEwtSet;\r
1439        //cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
1440 \r
1441        return(CV_SUCCESS);\r
1442      }\r
1443 \r
1444      private double N_VMin(NVector x) {\r
1445          int i;\r
1446          if ((x == null) || (x.data == null) || (x.data.length == 0)) {\r
1447                  return Double.NaN;\r
1448          }\r
1449          double data[] = x.data;\r
1450          int n = data.length;\r
1451          double min = data[0];\r
1452          for (i = 1; i < n; i++) {\r
1453                  if (data[i] < min) {\r
1454                          min = data[i];\r
1455                  }\r
1456          }\r
1457          return min;\r
1458      }\r
1459      \r
1460      /*\r
1461       * CVodeRootInit\r
1462       *\r
1463       * CVodeRootInit initializes a rootfinding problem to be solved\r
1464       * during the integration of the ODE system.  It loads the root\r
1465       * function pointer and the number of root functions, and allocates\r
1466       * workspace memory.  The return value is CV_SUCCESS = 0 if no errors\r
1467       * occurred, or a negative value otherwise.\r
1468       */\r
1469 \r
1470      int CVodeRootInit(CVodeMemRec cv_mem, int nrtfn, int g)\r
1471      {\r
1472        int i, nrt;\r
1473 \r
1474        /* Check cvode_mem */\r
1475        if (cv_mem==null) {\r
1476          cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeRootInit",\r
1477                         MSGCV_NO_MEM);\r
1478          return(CV_MEM_NULL);\r
1479        }\r
1480 \r
1481        nrt = (nrtfn < 0) ? 0 : nrtfn;\r
1482 \r
1483        /* If rerunning CVodeRootInit() with a different number of root\r
1484           functions (changing number of gfun components), then free\r
1485           currently held memory resources */\r
1486        if ((nrt != cv_mem.cv_nrtfn) && (cv_mem.cv_nrtfn > 0)) {\r
1487          cv_mem.cv_glo = null;\r
1488          cv_mem.cv_ghi = null;\r
1489          cv_mem.cv_grout = null;\r
1490          cv_mem.cv_iroots = null;\r
1491          cv_mem.cv_rootdir = null;\r
1492          cv_mem.cv_gactive = null;\r
1493 \r
1494          cv_mem.cv_lrw -= 3 * (cv_mem.cv_nrtfn);\r
1495          cv_mem.cv_liw -= 3 * (cv_mem.cv_nrtfn);\r
1496        }\r
1497 \r
1498        /* If CVodeRootInit() was called with nrtfn == 0, then set cv_nrtfn to\r
1499           zero and cv_gfun to NULL before returning */\r
1500        if (nrt == 0) {\r
1501          cv_mem.cv_nrtfn = nrt;\r
1502          cv_mem.cv_gfun = -1;\r
1503          return(CV_SUCCESS);\r
1504        }\r
1505 \r
1506        /* If rerunning CVodeRootInit() with the same number of root functions\r
1507           (not changing number of gfun components), then check if the root\r
1508           function argument has changed */\r
1509        /* If g != NULL then return as currently reserved memory resources\r
1510           will suffice */\r
1511        if (nrt == cv_mem.cv_nrtfn) {\r
1512          if (g != cv_mem.cv_gfun) {\r
1513            if (g < 0) {\r
1514              cv_mem.cv_glo = null;\r
1515              cv_mem.cv_ghi = null;\r
1516              cv_mem.cv_grout = null;\r
1517              cv_mem.cv_iroots = null;\r
1518              cv_mem.cv_rootdir = null;\r
1519              cv_mem.cv_gactive = null;\r
1520 \r
1521              cv_mem.cv_lrw -= 3*nrt;\r
1522              cv_mem.cv_liw -= 3*nrt;\r
1523 \r
1524              cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
1525                             MSGCV_NULL_G);\r
1526              return(CV_ILL_INPUT);\r
1527            }\r
1528            else {\r
1529              cv_mem.cv_gfun = g;\r
1530              return(CV_SUCCESS);\r
1531            }\r
1532          }\r
1533          else return(CV_SUCCESS);\r
1534        }\r
1535 \r
1536        /* Set variable values in CVode memory block */\r
1537        cv_mem.cv_nrtfn = nrt;\r
1538        if (g < 0) {\r
1539          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
1540                         MSGCV_NULL_G);\r
1541          return(CV_ILL_INPUT);\r
1542        }\r
1543        else cv_mem.cv_gfun = g;\r
1544 \r
1545        /* Allocate necessary memory and return */\r
1546        cv_mem.cv_glo = null;\r
1547        cv_mem.cv_glo = new double[nrt];\r
1548        if (cv_mem.cv_glo == null) {\r
1549          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1550                         MSGCV_MEM_FAIL);\r
1551          return(CV_MEM_FAIL);\r
1552        }\r
1553          \r
1554        cv_mem.cv_ghi = null;\r
1555        cv_mem.cv_ghi = new double[nrt];\r
1556        if (cv_mem.cv_ghi == null) {\r
1557          cv_mem.cv_glo = null;\r
1558          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1559                         MSGCV_MEM_FAIL);\r
1560          return(CV_MEM_FAIL);\r
1561        }\r
1562          \r
1563        cv_mem.cv_grout = null;\r
1564        cv_mem.cv_grout = new double[nrt];\r
1565        if (cv_mem.cv_grout == null) {\r
1566          cv_mem.cv_glo = null;\r
1567          cv_mem.cv_ghi = null;\r
1568          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1569                         MSGCV_MEM_FAIL);\r
1570          return(CV_MEM_FAIL);\r
1571        }\r
1572 \r
1573        cv_mem.cv_iroots = null;\r
1574        cv_mem.cv_iroots = new int[nrt];\r
1575        if (cv_mem.cv_iroots == null) {\r
1576          cv_mem.cv_glo = null;\r
1577          cv_mem.cv_ghi = null;\r
1578          cv_mem.cv_grout = null;\r
1579          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1580                         MSGCV_MEM_FAIL);\r
1581          return(CV_MEM_FAIL);\r
1582        }\r
1583 \r
1584        cv_mem.cv_rootdir = null;\r
1585        cv_mem.cv_rootdir = new int[nrt];\r
1586        if (cv_mem.cv_rootdir == null) {\r
1587          cv_mem.cv_glo = null; \r
1588          cv_mem.cv_ghi = null;\r
1589          cv_mem.cv_grout = null;\r
1590          cv_mem.cv_iroots = null;\r
1591          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1592                         MSGCV_MEM_FAIL);\r
1593          return(CV_MEM_FAIL);\r
1594        }\r
1595 \r
1596 \r
1597        cv_mem.cv_gactive = null;\r
1598        cv_mem.cv_gactive = new boolean[nrt];\r
1599        if (cv_mem.cv_gactive == null) {\r
1600          cv_mem.cv_glo = null; \r
1601          cv_mem.cv_ghi = null;\r
1602          cv_mem.cv_grout = null;\r
1603          cv_mem.cv_iroots = null;\r
1604          cv_mem.cv_rootdir = null;\r
1605          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1606                         MSGCV_MEM_FAIL);\r
1607          return(CV_MEM_FAIL);\r
1608        }\r
1609 \r
1610 \r
1611        /* Set default values for rootdir (both directions) */\r
1612        for(i=0; i<nrt; i++) cv_mem.cv_rootdir[i] = 0;\r
1613 \r
1614        /* Set default values for gactive (all active) */\r
1615        for(i=0; i<nrt; i++) cv_mem.cv_gactive[i] = true;\r
1616 \r
1617        cv_mem.cv_lrw += 3*nrt;\r
1618        cv_mem.cv_liw += 3*nrt;\r
1619 \r
1620        return(CV_SUCCESS);\r
1621      }\r
1622      \r
1623      class SUNLinearSolver {\r
1624          long N; // Size of the linear system, the number of matrix rows\r
1625          long pivots[]; // Array of size N, index array for partial pivoting in LU factorization\r
1626          long last_flag; // last error return flag from internal setup/solve\r
1627      }\r
1628      \r
1629      private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
1630          int matrixRows = A.length;\r
1631          int matrixColumns = A[0].length;\r
1632          if (matrixRows != matrixColumns) {\r
1633                  MipavUtil.displayError("Did not have the matrix rows = matrix columns required for a LinearSolver");\r
1634                  return null;\r
1635          }\r
1636          int vecLength = y.getData().length;\r
1637          if (matrixRows != vecLength) {\r
1638                  MipavUtil.displayError("Did not have the matrix rows equal to vector length required for a LinearSolver");\r
1639                  return null;\r
1640          }\r
1641          SUNLinearSolver S = new SUNLinearSolver();\r
1642          S.N = matrixRows;\r
1643          S.last_flag = 0;\r
1644          S.pivots = new long[matrixRows];\r
1645          return S;\r
1646      }\r
1647 \r
1648 \r
1649 }